Regularization Methods for Differential

Many mathematical models arising in science and engineering, including circuit and device simulation in VLSI, constrained mechanical systems in robotics and vehicle simulation, certain models in early vision and incompressible uid ow, lead to computationally challenging problems of di erential equations with constraints, and more particularly to high-index, semi-explicit di erential-algebraic equations (DAEs). The direct discretization of such models in order to solve them numerically is typically fraught with di culties. We thus need to reformulate the original problem into a better behaved problem before discretization. Index reduction with stabilization is one class of reformulations in the numerical solution of high index DAEs. Another class of reformulations is called regularization. The idea is to replace a DAE by a better behaved nearby system. This method reduces the size of the problem and avoids the derivatives of the algebraic constraints associated with the DAE. It is more suitable for problems with some sort of singularities in which the constraint Jacobian does not have full rank. Unfortunately, this method often results in very sti systems, which accounts for its lack of popularity in practice. In this thesis we develop a method which overcomes this di culty through a combination of stabilization and regularization in an iterative procedure. We call it the sequential regularization method (SRM). Several variants of the SRM which work e ectively for various circumstances are also developed. The SRM keeps the bene ts of regularization methods and avoids the need for using a sti solver for the regularized problem. Thus the method is an important improvement over usual regularization methods and can lead to improved numerical methods requiring only solutions to linear systems. The SRM also provides cheaper and more e cient methods than the ii usual stabilization methods for some choices of parameters and stabilization matrix. We propose the method rst for linear index-2 DAEs. Then we extend the idea to nonlinear index-2 and index-3 problems. This is especially useful in applications such as constrained multibody systems which are of index-3. Theoretical analysis and numerical experiments show that the method is useful and e cient for problems with or without singularities. While a signi cant body of knowledge about the theory and numerical methods for DAEs has been accumulated, almost none of it has been extended to partial di erential-algebraic equations (PDAEs). As a rst attempt we provide a comparative study between stabilization and regularization (or pseudo-compressibility) methods for DAEs and PDAEs, using the incompressible Navier-Stokes equations as an instance of PDAEs. Compared with stabilization methods, we nd that regularization methods can avoid imposing an arti cial boundary condition for the pressure. This is a feature for PDAEs not shared with DAEs. Then we generalize the SRM to the nonstationary incompressible Navier-Stokes equations. Convergence is proved. Again nonsti time discretization can be applied to the SRM iterations. Other interesting properties associated with discretization are discussed and demonstrated. The SRM idea is also applied to the problem of miscible displacement in porous media in reservoir simulation, speci cally to the pressure-velocity equation. Advantages over mixed nite element methods are discussed. Error estimates are obtained and numerical experiments are presented. Finally we discuss the numerical solution of several singular perturbation problems which come from many applied areas and regularized problems. The problems we consider are nonlinear turning point problems, a linear elliptic turning point problem and a second-order hyperbolic problem. Some uniformly convergent schemes with respect to the perturbation parameter are constructed and proved. A spurious solution phenomenon for the upwinding scheme is analyzed. iii Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgement ix 1 Introduction 1 1.1 Regularization for Di erential-Algebraic Equations (DAEs) . . . . . . 1 1.2 Regularization for the Incompressible Navier{ Stokes equations . . . . 7 1.3 A Problem in Reservoir Simulation . . . . . . . . . . . . . . . . . . . 8 1.4 Regularization for Di erential Equations without Constraints . . . . . 10 1.5 Contribution of This Thesis . . . . . . . . . . . . . . . . . . . . . . . 12 2 Sequential Regularization Methods for Di erential Algebraic Equations 15 2.1 Motivation of the SRM for General High Index DAEs . . . . . . . . . 15 2.2 Linear Index-2 Problems . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Problem Conditioning . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 Derivation of the SRM . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.5 Convergence Analysis of the SRM . . . . . . . . . . . . . . . . . . . . 30 2.6 Discretization and Implementation Issues . . . . . . . . . . . . . . . . 34 2.7 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.8 More about the Proof of Theorem 2.1 . . . . . . . . . . . . . . . . . . 46 iv 3 SRM for Nonlinear Problems 52 3.1 Nonlinear, Nonsingular Index-2 Problems . . . . . . . . . . . . . . . . 53 3.1.1 The case 1 = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.1.2 The case 1 = 0 . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.2 Nonlinear, Singular Index-2 Problems . . . . . . . . . . . . . . . . . . 62 3.3 SRM for Nonlinear Higher-index Problems . . . . . . . . . . . . . . . 64 3.3.1 The case of nonsingular GB . . . . . . . . . . . . . . . . . . . 65 3.3.2 The case for constraint singularities . . . . . . . . . . . . . . 71 3.4 SRM for Constrained Multibody Systems . . . . . . . . . . . . . . . . 72 3.5 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4 SRM for the Nonstationary Incompressible Navier-Stokes Equations 834.1 DAE Methods for Navier-Stokes Equations . . . . . . . . . . . . . . . 83 4.2 Preliminaries and the Properties of the Regularized Problems . . . . 88 4.3 Convergence of the SRM . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.3.1 Two linear auxiliary problems . . . . . . . . . . . . . . . . . . 94 4.3.2 The error estimate of SRM . . . . . . . . . . . . . . . . . . . . 97 4.4 Discretization Issues and Numerical Experiments . . . . . . . . . . . 103 5 SRM for the Simulation of Miscible Displacement in Porous Media112 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.2 SRM Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.3 Convergence Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.4 The Galerkin Approximation and Its Error Estimates . . . . . . . . . 120 5.5 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6 Numerical Methods of Some Singular Perturbation Problems 127 v 6.1 One Dimensional Quasilinear Turning Point Problems . . . . . . . . . 128 6.1.1 A repulsive turning point problem . . . . . . . . . . . . . . . . 128 6.1.2 An attractive turning point problem . . . . . . . . . . . . . . 133 6.2 Notes about Spurious Solutions of Upwind Schemes . . . . . . . . . . 139 6.2.1 Inadequacy of Yavneh's argument . . . . . . . . . . . . . . . . 139 6.2.2 Our explanation . . . . . . . . . . . . . . . . . . . . . . . . . . 141 6.3 A Linear Hyperbolic-Hyperbolic Singularly Perturbed Initial-Boundary Value Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 6.3.1 Construction of asymptotic solution and its remainder estimate 147 6.3.2 Di erence scheme and its uniform convergence . . . . . . . . . 150 7 Conclusion and Future Work 154 7.1 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . 154 7.2 Discussion of future work . . . . . . . . . . . . . . . . . . . . . . . . . 157 Bibliography 160 vi List of Tables 2.1 SRM errors for Example 2.2 using the midpoint scheme . . . . . . . 43 2.2 SRM errors for Example 2.2 using the shooting-back technique . . . 44 2.3 SRM errors for Example 2.3 using backward Euler . . . . . . . . . . 44 2.4 SRM errors for Example 2.3 using forward Euler . . . . . . . . . . . 45 2.5 Errors near singularity using modi ed formula (2.42) . . . . . . . . . 45 2.6 Errors for problem without singularity using modi ed formula (2.42) 45 3.1 Errors for Example 3.3 using the explicit second order Runge-Kutta scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.2 Example 3.4 { bounded y and singularity at t = :5 . . . . . . . . . . 78 3.3 Example 3.5 { unbounded y and singularity at t = :5 . . . . . . . . 78 3.4 Errors for Example 3.6 using SRM (3.45)-(3.46) . . . . . . . . . . . . 80 3.5 Drifts of the SRM for the slider-crank problem . . . . . . . . . . . . 81 4.1 SRM errors for = 0:1 without upwinding . . . . . . . . . . . . . . . 109 4.2 SRM errors for = 0:001 with upwinding . . . . . . . . . . . . . . . 110 4.3 SRM errors for = 0:001 with a pretty large time step k = h = 0:1 . 110 4.4 SRM errors for = 0:1 with 1 = 0 . . . . . . . . . . . . . . . . . . 111 5.1 Numerical results for Example 5.1 with grid size = 1 40 . . . . . . . . 126 5.2 Numerical results for Example 5.1 with grid size = 1 40 . . . . . . . . 126 5.3 Numerical results for Example 5.2 . . . . . . . . . . . . . . . . . . . 126 vii List of Figures 2.1 planar slider-crank: initial state in solid line, subsequent states in dotted lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1 Two-link planar robotic system . . . . . . . . . . . . . . . . . . . . . 79 3.2 Solution for slider-crank problem with singularities . . . . . . . . . . 81 3.3 Acceleration of slider end . . . . . . . . . . . . . . . . . . . . . . . . 82 5.1 One element with velocity on each edge and pressure at the center . 123 viii Acknowledgement I am specially grateful to Dr. Uri Ascher, my thesis supervisor, who contributed signi cantly to the development of this thesis. That contribution ranged from fruitful suggestions and enlightening discussions to general guidelines and ideas about how to present the results. I am also grateful to Drs. John Heywood, Tim Salcudean, Michael Ward and Brian Wetton, members of my supervisory committee, for their helpful comments and suggestions during this research and writing of the thesis. I would also like to thank the nancial supports from Killam pre-doctoral fellowship committee and fromDepartment of Mathematics, University of British Columbia. Most of all, I wish to acknowledge the tolerance displayed by my wife, Yuchun, and by my daughter, Xueshu (Susan) during this work. This thesis work often dominated the evenings, weekends and summers which could have been used in activities of greater interest to them. ix Chapter 1 Introduction Many mathematical models arising in science and engineering, including circuit and device simulation in VLSI, constrained mechanical systems in robotics and vehicle simulation, certain models in early vision and incompressible uid ow, lead to computationally challenging problems of di erential equations with constraints, and more particularly to high-index, semi-explicit di erential-algebraic equations (DAEs). The direct discretization of such models in order to solve them numerically is typically fraught with di culties, and most methods proposed in the literature seek to circumvent this by employing combinations of problem reformulation, regularization 1 and special discretization techniques. We will consider the regularization of such mathematical models and the numerical solution of the resulting regularized formulations. These formulations are often singular perturbation problems because they typically depend on a small parameter which provides a measure of the closeness between the regularized and the original problems. We will also apply our regularization method and idea to other relevant practical problems. 1.1 Regularization for Di erential-Algebraic Equations (DAEs) DAEs are special implicit ordinary di erential equations (ODEs) f(x0(t); x(t); t) = 0; (1.1) 1The concept of regularization was introduced by Tikhonov (see [109]). Its idea is that one solves a better behaved nearby problem instead of solving the original problem to circumvent some sort of di culties. See the next section for more details.1 Chapter 1. Introduction 2 where the partial Jacobian matrix fy(y; x; t) is singular for all relevant values of its arguments. Here x0 = dx dt . An extension to partial di erential equations is considered in the next subsection. DAEs were motivated by applications like network analysis, circuit simulation and mechanical system simulation starting in the 1970's. They often arise as ordinary di erential equations with additional variables and (equality) algebraic constraints. An extensive list of applications is given in [92]. In the 1980's, DAEs have developed into a highly topical subject of computational and applied mathematics. Contributions devoted to DAEs have appeared in various elds, such as applied mathematics, scienti c computation, mechanical engineering, chemical engineering, system theory, etc. Frequently, other names have been assigned to DAEs, e.g. semistate equations, descriptor systems, singular systems. Gear ([51], 1971) proposed to handle DAEs numerically by backward di erentiation formulas (BDF). For a long time DAEs had been considered not to di er essentially from regular implicit ODEs in general. Only since about 1980, because of computational results that could not be brought into line with the above supposition (e.g. Sincovec et al [103], 1981), the mathematical and particularly the numerical community have started investigating DAEs more thoroughly. With their famous paper, Gear, Hsu and Petzold ([52], 1981) started a discussion on DAEs that will surely be carried on for a while. The structure of DAEs is very much related to the concept of index, which is a measure of the amount of singularity of the system. There are several ways to de ne index. The most popular one is called di erential index, which is de ned as the minimal number of analytical di erentiations in t such that (1.1) can be transformed by algebraic manipulations into an explicit ordinary di erential system (in the original unknowns) x0 = (x; t) Chapter 1. Introduction 3 which in turn is called the \underlying ODE ". There are several structural forms of DAEs which appear frequently in applications (see [92]). The di erential index of these structural forms can be found by di erentiating their algebraic constraints with respect to t (and substituting into the di erential equations which complement the algebraic constraints). For instance, Semi-explicit index-1 system x0 = f(x; y); (1.2a) 0 = g(x; y); (1.2b) if gy is invertible. Hessenberg index-2 system x0 = f(x; y); (1.3a) 0 = g(x); (1.3b) if gxfy is invertible. Hessenberg index-3 system x0 = f(x; y); (1.4a) y0 = k(x; y; z); (1.4b) 0 = g(x); (1.4c) if gxfykz is invertible. Mechanical multibody systems with holonomic constraints are examples of Hessenberg index-3 DAEs. In [56] it is pointed out that higher index ( 2) DAEs, in the natural function space formulations, lead to ill-posed 2 problems because they do not have the usual 2A problem is called ill-posed if it is not well-posed. A problem is called well-posed if it satis es three conditions, i.e. the existence, uniqueness and stability of its solution. \Stability" means that the solution of the problem continuously depends on the \data", which may be initial data, boundary data, coe cients in the equation, values of the operator, etc.. Here high-index DAEs fail to satisfy a stability condition. Chapter 1. Introduction 4 stability property of di erential equations in general. Example 1.1 Consider (See [86])x0 = y; (1.5a) 0 = x p(t); (1.5b) where x(t) and y(t) are scalar functions and p(t) is a given function. This is a very simple index-2 DAE. The exact solution is x = p(t), y = p0(t). If we add a small perturbation sin!t, << 1, to the right hand side of the second equation we have the exact solution x = p(t) + sin!t; y = p0(t) + ! cos!t: Hence y y = ! cos!t could be very large if ! 1= , i.e. the solution changes a lot under a small change in the right-hand side of the equation. 2 A numerical method which is directly applied to a complex, ill-posed problem may generally fail. Therefore, to solve DAEs, we have to stabilize such problems to bring about continuous dependence on the \data" (or stability). One such approach is to change the formulation of the problem but not its solution, e.g. in Example 1.1 we di erentiate (1.5b) once and obtain x0 = y (1.6a) 0 = y p0(t); x(0) = p(0): (1.6b) Now (1.6) is of index-1 and becomes a well-posed problem with the same solution as (1.5). We can solve (1.6) instead of (1.5) and gain well-posedness. However, such a direct index reduction procedure may cause the well-known drift di culty (see [29]), i.e. the approximate solution of (1.6) may be far from satisfying the constraint (1.5b). Chapter 1. Introduction 5 Hence, methods have been designed to prevent moving away from the constraints. Baumgarte's stabilization [17] and projection invariant methods (cf. [16]) are popular among such methods. Most of these approaches treat initial value problems, and only a few apply to boundary value problems. See [29, 58] for various numerical methods for initial value problems; [93, 7, 15] for boundary value problems; and [16, 9, 35] for a survey of various stabilization reformulations. Another approach consists in adding some small perturbation terms (measured by a small positive parameter ) to the given DAE. The perturbed problem is close to the original problem (if is small) and is well-posed. Such an approach is usually called regularization. This is a natural approach since the high-index DAE is illposed; indeed in [43] a well-known Tikhonov regularization algorithm (see [109]) was applied to solve DAEs. However, such a method seems so general that it is not su ciently related to the special structure of the DAE. There are two types of regularization methods which are probably more interesting for DAE researchers. One is called parameterization. One such possibility, the pencil regularization, was given independently by Boyarintsev [24] (or see his newer book [25] published in English) and Campbell [32]. But the regularized problem is ensured to be well-posed only for constant coe cient cases. A further parameterization was proposed by M arz [85]. Her regularization is aimed at obtaining well-posed index-1 DAEs instead of obtaining well-posed ODEs. Heuristically, it seems evident that the DAE is less changed if it is transformed into an index-1 DAE rather than an ODE. Marz's regularization was proved to be well-posed for usual structural forms of DAEs. We refer to [59] for further results in this direction. Another class of regularization uses the penalty idea (see [84, 91]). It originates from penalty methods for constrained optimization problems. Note that an algebraic equation in a DAE can be viewed as a constraint. This method seems more natural for DAEs. References [68, 70, 69] used the penalty regularization and singular Chapter 1. Introduction 6 perturbation theory to determine the solutions of DAEs when the initial or boundary values are given improperly (i.e. inconsistently). In Chapter 2 we will mention these methods again with a bit more detail and indicate that Marz's regularization is actually a kind of penalty method. Because the regularization method requires fewer di erentiations of the constraints it is perhaps more suitable for DAEs which have singularities, i.e. whose constraints do not have full rank, e.g. when the matrix gxfy is singular at some isolated points in the index-2 system (1.3). These problems can be challenging for the methods that are usually employed and appear frequently in simulation of constrained mechanical systems. To our knowledge, there has not been a paper in the numerical analysis literature about this until the recent two preprints [11] and [94] (although a number of relevant papers appear in the mechanical engineering literature). From a practical point of view, a number of codes which work well and e ciently (at least if the regularization parameter is not too small) are available for numerically solving the regularized problems. We also note that the regularization method requires less smoothness of the coe cients of the di erential-algebraic problem than other stabilization methods. These are the advantages of the regularization method. The dominant disadvantage in the above regularization methods is that the parameter must be small enough to maintain the accuracy of the numerical method we use for the regularized problem at an acceptable level. Hence, a sti solver is necessary. Typically in regularization methods, the parameter must be chosen both \large enough" and \small enough": large so that the regularized problem would behave signi cantly better than the original, and small so that its solution will not di er too much from that of the original problem. We will present a class of new regularization methods, inspired by [19], which we call sequential regularization method (SRM). The SRM can be viewed as a combination of the penalty method with Baumgarte's stabilization in an iterative procedure; Chapter 1. Introduction 7 see x2.4 or x3.1 for speci c instances. It is applicable for DAEs with constraint singularities. Moreover, the regularization parameter in the method is not necessarily small. Thus, a nonsti solver can be used for solving the regularized problems. Some variants of the SRM are discussed for index-2 and index-3 DAEs with the goal of simplifying the computations. We will apply the SRM to mechanical multibody systems as well. 1.2 Regularization for the Incompressible Navier{ Stokes equations As noted before, DAEs have become a highly topical subject of applied and numerical mathematics. However, there seems to be still a void in the literature about partial di erential equations with constraints (PDAEs). A typical instance of such problems is the well-known incompressible Navier-Stokes equations: ut + (u grad)u = u gradp+ f ; (1.7a) divu = 0; (1.7b) uj@ = b ; ujt=0 = a; (1.7c) in a bounded twoor three-dimensional domain and 0 t T . Here u(x; t) represents the velocity of a viscous incompressible uid, p(x; t) the pressure, f the prescribed external force, a(x) the prescribed initial velocity, and b(t) the prescribed boundary values. The system (1.7) can be seen as a partial di erential equation with constraint (1.7b) with respect to the time variable t. Hence, we call it a PDAE. It is easily veri ed that it is of index-2 without singularities since the operator divgrad = is invertible (under appropriate boundary conditions). A huge number of methods have been designed to solve the nonstationary incompressible Navier-Stokes equations (1.7). Direct discretizations include nite di erence Chapter 1. Introduction 8 and nite volume techniques on staggered grids (e.g. [65, 26, 66]), nite element methods using conformal and nonconformal elements (e.g. [54, 110, 63, 64]) and spectral methods (e.g. [33]). Another approach yielding many methods has involved some initial reformulation and/or regularization of the equations, to be followed by a discretization of the (hopefully) simpli ed system of equations. Examples of such methods include pseudo-compressibility methods, projection and pressure-Poisson reformulations (e.g. [36, 55, 72, 97, 102, 117]). The two types of regularizations we mentioned in x1.1 for DAEs were already proposed in the Navier-Stokes context quite a while ago (cf. [108, 72]). We are interested in the generalization of the SRM to this problem because the regularized problems can be made essentially nonsti and then a more convenient di erence scheme (e.g. an explicit scheme) in time is possible. Moreover, the method retains the bene ts of the penalty method. For example, computations for the velocity u and the pressure p are uncoupled and an arti cial boundary condition for calculating the pressure p is not necessary. 1.3 A Problem in Reservoir Simulation The idea of the SRM can be applied to a reservoir-simulation problem | miscible displacement in porous media. Miscible displacement is an enhanced oil-recovery process that has attracted considerable attention in the petroleum industry. It involves injection of a solvent at certain wells in a petroleum reservoir, with the intention of displacing resident oil to other wells for production. This oil may have been left behind after primary production by reservoir pressure and secondary production by water ooding. The economics of the process can be precarious, because the chemicals it requires are expensive and the performance of the displacement is by no means guaranteed. Complex physical behavior in the reservoir will determine whether enough additional oil is recovered Chapter 1. Introduction 9 to make the expense worthwhile. A numerical simulation of the complex process undoubtedly plays an important role. Mathematically, the process is described by a convection{ dominated parabolic partial di erential equation for each chemical component in the system. By summing up the component equations, one can obtain an equation that determines the pressure in the system; this nonlinear equation is elliptic or parabolic, according to whether the system is incompressible or compressible. Thus, in this problem one encounters elliptic, parabolic, and near{hyperbolic equations with complicated nonlinear behavior.For simplicity, we consider the miscible displacement of one incompressible uid by another in a porous reservoir R2 over a time period [0; T ]. Let p(x; t) and u(x; t) denote the pressure and Darcy velocity of the uid mixture, and c(x; t) the concentration of the invading uid. Then the mathematical model is a strongly coupled nonlinear system of partial di erential equations (see [44, 82]): u = a(gradp gradd); (x; t) 2 [0; T ]; (1.8a) divu = q(x; t); (x; t) 2 [0; T ]; (1.8b) @c @t div(D(u)gradc) + u gradc = g(c); (x; t) 2 [0; T ]; (1.8c) with the boundary conditions u n = 0; (x; t) 2 [0; T ]; (1.9a) D(u)gradc n = 0; (x; t) 2 [0; T ]; (1.9b) and initial condition c(x; 0) = c0(x); x 2 ; (1.10) where a = a(x; c) is the mobility of the uid mixture, = (x; c) and d(x) are the gravity and vertical coordinate, q is the imposed external rates of ow, (x) is Chapter 1. Introduction 10 the porosity of the rock, D is the coe cient of molecular di usion and mechanical dispersion of one uid into the other, g = g(x; t; c) is a known linear function of c representing sources, and n is the exterior normal to the boundary = @ . The pressure-velocity equation (1.8a){(1.8b) is elliptic (after eliminating u). The concentration equation (1.8c) is parabolic, but normally convectiondominated. It is derived from the conservation of mass which involves the Darcy velocity of the uid mixture, but the pressure variable does not appear in it. Thus a good approximation of the concentration equation requires accurate solution for the velocity variables. Mixed nite element methods have been applied to the pressure-velocity equation, which can yield velocity solutions one order more accurate than those obtained using corresponding nite di erence and usual nite element methods [40, 41, 42, 46, 47, 120]. However, the nite dimensional spaces for the velocity and pressure need to satisfy the Babuska-Brezzi condition, and the resulting linear system does not have a positive de nite coe cient matrix. Moreover, the number of degrees of freedom in the linear system doubles that of nite di erence or nite element methods. We are interested in designing an SRM for the pressure-velocity equation since the SRM formulation can produce as accurate a velocity approximation as the mixed nite element methods, and can avoid the above-mentioned problems in mixed nite element methods. 1.4 Regularization for Di erential Equations without Constraints Regularization methods are also used to treat di erential equations without constraints when these equations have some sort of singularities or their solutions may have some sort of discontinuities. Examples are viscous solutions of hyperbolic conservation laws [76], shape from shading problems with singularities [34] and transition Chapter 1. Introduction 11 phenomena in semi-conductor device simulation [13]. A frequently considered problem is the following rst-order partial di erential equation a(u; x; t)@u @t + n Xi=1 bi(u; x; t) @u @xi + c(u; x; t) = 0 (1.11) with appropriate initial and boundary conditions. In general, (1.11) may have some kind of discontinuous solutions, e.g. a shock wave [76]. Or, if we consider the steadystate case, (1.11) may have singularities, i.e. bi may vanish at some points, as in the shape from shading problem [34]. Some regularization techniques have been designed for solving (1.11). A popular one is a(u; x; t)@u @t + d(u; x; t) u+ n Xi=1 bi(u; x; t) @u @xi + c(u; x; t) = 0: (1.12) The regularized problem often has a physical meaning by itself, e.g. a time-dependent advection-di usion equation with a small di usion term. Another choice could be d(u; x; t)(@2u @t2 u) + a(u; x; t)@u @t + n Xi=1 bi(u; x; t) @u @xi + c(u; x; t) = 0: (1.13) This has physical meanings as well, e.g. a tra c ow problem [118] or so-called overdamped vibration problems [104]. Thus the approximate resolution of (1.11) becomes that of the singular perturbation problem (1.12) or (1.13). Unlike the SRM, the regularization parameter has to be small in comparison with the mesh size to ensure that the solution of the regularized problem be a good approximation of that of the original problem. It is well-known that there are di culties in solving these regularized problems numerically with small , e.g. the stability problem for the central di erence scheme and the accuracy problem for the upwinding scheme in a boundary layer region in which the derivatives of the solution may be large (see [37, 61]). We will consider some special cases of (1.12) or (1.13) and focus mainly on uniformly convergent methods. We will discuss spurious solutions of a simple upwinding scheme as well. Chapter 1. Introduction 12 1.5 Contribution of This Thesis Our objectives are to propose and to investigate regularization methods for various di erential equations with or without constraints. Most attention is paid to ordinary and partial di erential equations with constraints (DAEs and PDAEs). We propose and analyze a regularization method called the sequential regularization method (SRM). A very important advantage of our regularization method (SRM) is that the problem after regularization need not be sti . Hence explicit di erence schemes can be used to avoid solving nonlinear systems and they make the computation much simpler. Improvements over stabilization methods and extra bene ts for PDAEs are also achieved. In Chapter 2, the SRM is presented for linear index-2 DAEs with or without constraint singularities. A complete theoretical analysis is performed for both cases. It is proved that the di erence between the exact solution of a DAE and the corresponding iterate becomes O( s) in magnitude at the sth iteration, at least away from the starting value of the independent variable t. Hence, the regularization parameter need not be very small so the regularized problems are less sti . By some choice of parameters the regularized problems can be essentially nonsti for any . As an example, a simple di erence scheme for solving the regularized problems is investigated. Implementation techniques are discussed to get an approximation in the whole region for boundary value problems and to economize storage for initial value problems. Numerical experiments support our theoretical results. Numerical examples also show that usual stabilization methods do not work for problems with constraint singularities. Most parts of this chapter are taken from the paper [11]. In Chapter 3, we extend the SRM to nonlinear problems and to DAEs with index higher than 2. Again, nonsti ness of the regularized problems is achieved. Rather Chapter 1. Introduction 13 than having one \winning" method, this is a class of methods from which a number of variants are singled out as being particularly e ective methods in certain circumstances of practical interest. In the case of no constraint singularity we prove convergence results. The method is also applied to constrained multibody systems. Numerical experiments con rm our theoretical predictions and demonstrate the viability of the proposed methods. Most parts of this chapter are taken from the paper [12]. In Chapter 4, we generalize the SRM to PDAEs, in particular, to the nonstationary Navier-Stokes equations. The convergence rate O( s) at the sth iteration is again proved for this PDAE case in appropriate norms. The SRM not only avoids the sti ness of the regularized problems which always occurs in pseudo-compressibility methods but also avoids providing an unphysical pressure boundary condition which has to be imposed in stabilization methods. Discretization and implementation issues of the SRM are considered as well. In particular, a simple explicit di erence scheme is analyzed and its stability is proved under the usual step size condition (independent of the regularization parameter ) of explicit schemes. The stability result also indicates that the step size restriction can be relaxed as the viscosity becomes small. A numerical example is calculated to demonstrate these results. The SRM formulation is new in the Navier{Stokes context and it performs well. Most parts of the chapter are taken from the paper [79]. In Chapter 5, we apply the idea of the SRM to the simulation of miscible displacement in porous media. The problem is modeled by a nonlinear coupled system of two partial di erential equations: the pressure-velocity equation and the concentration equation. Only the approximation of velocity is important for the approximation of concentration. The SRM idea is used for the pressure-velocity equation. An O( s) error estimate at the sth SRM iteration is also proved. A Galerkin nite element Chapter 1. Introduction 14 method is used for the discretization of the SRM formulation. It is capable of producing as accurate a velocity approximation as the mixed nite element method. But unlike the mixed nite element method its sti ness matrix is symmetric positive de nite and its nite element spaces need not satisfy the so-called Babuska-Brezzi condition. Most parts of this chapter are taken from the report [82]. Chapter 6 is devoted to numerical methods of some special cases of singular perturbation equations in the form of (1.12) or (1.13). Sections x6.1 and x6.3 describe a collection of papers [79, 115, 107] which re ects the author's earlier research interests. We believe that, by considering these special cases, we make steps towards the general problems (1.12) or (1.13) which are undoubtedly very di cult. Also, these special cases have practical meaning in themselves, hence, are worthwhile to be considered independently. In x6.1, we consider the one dimensional steady-state case of (1.12). x6.2 covers a special two-dimensional steady-state instance of (1.12) given in [121] to show that upwind schemes may lead to spurious solutions even for problems with very smooth solutions. We indicate that this is actually an ill-posed problem when is small. Hence, it is not strange that a direct discretization to the problem fails. In x6.3, we consider the linear one dimensional time-dependent case of (1.13). In this case, derivatives of the reduced problemmay be discontinuous along the characteristic curves. Finally, conclusions and possible future work are contained in Chapter 7. Chapter 2 Sequential Regularization Methods for Di erential Algebraic Equations 2.1 Motivation of the SRM for General High Index DAEs The sequential regularization method (SRM) is motivated from the augmented Lagrangian method applied to constrained multibody systems (index-3 DAEs in general) by Bayo and Avello [19] and an earlier paper [20]. So we start this chapter by considering a mechanical system whose con guration is characterized by the generalized coordinates q. Let L be the system Lagrangian, de ned by L = T V; where T and V are the kinetic energy and the potential energy, respectively. Let Q represent non-conservative forces. Usually the Lagrangian coordinates are not independent, but rather are interrelated through certain constraint conditions. When the connections between bodies are of holonomic type, these constraint conditions can be expressed mathematically in the following form: (q; t) = 0: (2.1) Then Hamilton's principle leads to the Euler-Lagrange equations: d dt(@L @q0) @L @q + Tq = Q; (2.2) where is a vector function whose components are Lagrange multipliers. For general multibody systems, (2.2) becomes M(q)q00+ Tq = f(q; q0) (2.3) 15 Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations16 with common initial conditions, where M is the mass matrix and q is the Jacobian of the constraint equations. (2.3) and (2.1) form the Euler-Lagrange equations for a constrained multibody system. This is an index-3 DAE if q has full rank. We have indicated that direct discretization would not be good in general for such a higher index DAE. One usual way to treat this problem is index reduction (to an index-1 DAE or an ODE). The most straightforward transformation of the DAE (2.3), (2.1) to an index-1 DAE involves replacing the constraint (2.1) with its second derivative plus initial conditions: 00 = d2 (q(t); t) dt2 = 0; (2.4a) (q(0); 0) = d (q(0); 0) dt = 0; (2.4b) (cf. (1.6)). However, this causes well-known drift di culties, i.e. the numerical solution of (2.3), (2.4) may drift away from the original constraints (2.1) as time proceeds. Hence we have to look for stabilized index reduction methods. A very popular method called Baumgarte's method proposed in 1972 [17] is a generalization of (2.4). It replaces (2.4a) with the equation 00 + a 0 + b = 0; (2.5) where a and b are parameters chosen so that the roots of the polynomial ( ) = 2 + a + b = 0; (2.6) are both negative, i.e. the initial value problem for the di erential equation (2.5) for is asymptotically stable (see [16]). The system (2.5), (2.3) can be written in the form 24 M Tq q 0 3524 q00 35 = 24 f(q; q0) ( q)0q0 ( t)0 a 0 b 35 : (2.7) The matrix 24 M Tq q 0 35 (2.8) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations17 -1.5 -1 -0.5 0 0.5 1 1.5 0 1 2 Figure 2.1: planar slider-crank: initial state in solid line, subsequent states in dotted lines is nonsingular if q has full row rank. Hence (2.7) can be integrated using standard numerical integrators. If q is rank-de cient, we have a potential di culty in solving (2.7). Baumgarte's method may not work then. The problem is called singular if q does not have full rank. Example 2.1 Consider two linked bars (see Fig. 2.1). One end of one bar is xed at the origin, allowing only rotational motion in the plane. The other end of the other bar slides on the x-axis. The equations of motion form a nonlinear index-3 DAE p0 = v Mv0 = f GT g(p) = 0 where xi; yi; i are the coordinates of the center of mass of the ith bar, and Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations18 p = (x1; y1; 1; x2; y2; 2)T : If the left bar is strictly shorter than the right bar, then the Jacobian matrix G of the constraint functions of this problem has full rank. The problem is nonsingular. If the length of these two bars are the same, for example, each with length 2 and mass 1, then we have M = diagf1; 1; 1=3; 1; 1; 1=3g f = (0; 9:81; 0; 0; 9:81; 0)T g = 0BBBBBBBBBBB@ x1 cos 1 y1 sin 1 x2 2x1 cos 2 y2 2y1 sin 2 2 sin 1 + 2 sin 2 1CCCCCCCCCCCA G = gp = 0BBBBBBBBBBB@ 1 0 sin 1 0 0 0 0 1 cos 1 0 0 0 2 0 0 1 0 sin 2 0 2 0 0 1 cos 2 0 0 2 cos 1 0 0 2 cos 2 1CCCCCCCCCCCA Clearly, as the mechanism moves left through the point where both bars are upright ( 1 = 2 ; 2 = 3 2 ) the last row of G vanishes at this one point and a singularity is obtained. When arriving at this point with no momentum, this is actually a bifurcation point where two subsequent motion con gurations are possible. We will consider only the case where the sliding bar continues to slide along the x-axis past the singularity, and note that the solution is smooth in the passage through the singularity. 2 In [19], Bayo and Avello proposed to solve the multibody system (2.3) using an augmented Lagrangian algorithm which is transplanted from the same method in the optimization context [5]. Their idea is to derive a modi ed formulation by adding to the expression of Hamilton's principle three terms: Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations19 a ctitious potential V =Xk 12 k!2 k 2k (2.9) a set of Rayleigh dissipative forces Gk = 2 k!k k 0k (2.10) a ctitious kinematic energy term T =Xk 1 2 k 02 k ; (2.11) where each k is a very large real number (the penalty), and !k and k represent the natural frequency and the damping ratio of the penalized system (mass, dashpot and spring) corresponding to the constraint k = 0. Then we get a modi ed EulerLagrange equation Mq00 + Tq ( 00 + 2 0 + 2 ) + Tq = f(q; q0) (2.12) or(M + Tq q)q00+ Tq = f(q; q0) Tq (( Tq )0q0 + ( t)0 + 2 0 + 2 ); (2.13) where , and are diagonal matrices that contain the values of the penalties, the natural frequencies and the damping ratios of the ctitious penalty systems assigned to each constraint condition. Because is not given in advance Bayo, Jalon and Serna [20] propose an iteration to solve (2.12) or (2.13) by comparing (2.12) with (2.3): s = s 1 + ( 00 + 2 0 + 2 )jq=qs 1; s = 1; 2; (2.14) to get an approximation of . In [19] the authors called the whole procedure an augmented Lagrangian algorithm and claimed that the algorithm works well, however without any theoretical analysis. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations20 In fact, the algorithm would certainly not work in general. In multibody motions q usually remains smooth even in passage through singular positions. However, as indicated in later sections, may be unbounded at the singularity. So the iteration (2.14), as an approximate procedure to obtain , is not appropriate in some cases. Also, 00 should not be included in (2.12) in the singular case because unbounded coe cients may appear in it. From (2.14) the formulation corresponds to Baumgarte's formulation since when s gets close to s 1 (2.14) gets close to Baumgarte's stabilization. But Baumgarte's stabilization is not as good as some other stabilization techniques (see [8]). Moreover, in [19] the authors indicated that the iteration (2.14) is applied until kq00 s q00 s 1k < , where is a user-speci ed tolerance. This criterion is perhaps applied because they did not know the convergence rate of their iterative procedure. We do not recommend this criterion because it causes not only unnecessary extra iterations but also makes a storage-saving implementation di cult (cf. x2.6). Our aim is to construct a method for general DAEs which not only avoids the above shortcomings, but for which we also do not have a Lagrangian and Hamilton's principle. So another derivation of the algorithm is needed. In this and the next chapter, we propose a class of algorithms motivated by the augmented Lagrangian method for more general DAEs of order , x( ) = f(x; x0; : : : ; x( 1); t) B(x; t)y; (2.15a) 0 = g(x; t): (2.15b) The DAE (2.15) has index +1 if GB is nonsingular for all t, 0 t tf , where G = gx. We are interested in the cases = 1 or 2. The Euler-Lagrange equations (2.3) for mechanical systems with holonomic constraints are in this form with = 2. The algorithm is derived by combining a modi ed penalty idea (a kind of regularization) given in [91] with stabilization techniques such as Baumgarte's stabilization or the stabilization analyzed in [8, 35] in an iterative procedure. We call the method the Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations21 sequential regularization method (SRM) (cf. [11, 12]). The method is applicable for more general higher index DAEs. More importantly, it works for both boundary and initial value problems and is justi ed by a theoretical analysis. The number of iterations can be determined beforehand depending on the penalty parameter = 1= , the mesh size h and the order of the method used. Since we specify the iteration number in advance we can design a procedure for the initial value case to perform the iteration without the need to store all previous approximate values. The sequential regularization method is actually a functional iteration procedure in which the di erence between the exact solution of a DAE and the corresponding iterate becomes O( s) in magnitude at the sth iteration, at least away from the starting value of the independent variable (which we shall call `time'). Hence, unlike the usual regularization, the perturbation parameter does not have to be chosen very small, so the regularized problems can be less sti and/or more stable. Next we will propose and analyze the SRM for the linear index-2 case with singularities. Numerical experiments are given to verify our theoretical results. Some simple di erence schemes for the regularized problems and implementation issues for the SRM are also discussed. 2.2 Linear Index-2 Problems We rst write down the linear index-2 problem: x0 = A(t)x B(t)y + q(t); (2.16a) 0 = G(t)x+ r(t) g(x; t); (2.16b) where A(t), B(t) and G(t) are su ciently smooth functions of t; 0 t tf , A(t) 2 Rnx nx , B(t) 2 Rnx ny , G(t) 2 Rny nx , and ny nx. We consider the DAE (2.16) subject to nx ny boundary conditions Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations22 B0x(0) + B1x(tf) = : (2.17) These boundary conditions are assumed to be such that they yield a unique solution for the ODE (2.16a) on the manifold given by (2.16b). In particular, if we were to replace (2.16b) by its di erentiated form 0 = Gx0 +G0x+ r0 = d dtg(x; t); (2.18a) g(x(0); 0) = G(0)x(0) + r(0) = 0; (2.18b) and use (2.18a) in (2.16a) to eliminate y and obtain nx ODEs for x, then the boundary value problem for x with (2.17) and (2.18b) speci ed has a unique solution. In the initial value case B1 = 0, this means that (2.17) and (2.18b) can be solved uniquely for x(0). We will give a more precise assumption in Lemma 2.1 below. The problem (2.16), (2.17) is of index-2 if GB is nonsingular for all t. However, here we allow GB to be singular at some isolated points of t. For simplicity of exposition, let us say that there is one singular point t ; 0 < t < tf . The inhomogeneities are q(t) 2 Rnx and r(t) 2 Rny . We are only interested in the kind of singularities as in Example 2.1, where the solution x of (2.16), (2.17) passes through the singularity smoothly. Returning to Example 2.1 (where B =M 1GT ), we can verify that, although the matrix GM 1GT is singular at the singularity , the matrix M 1GT (GM 1GT ) 1G is smooth for all t. Also, two types of singular constraints (i.e., with vanishing rows or with some rows linearly dependent at some points) mentioned in [2] both have a similar property. Thus, for the linear model (2.16), we assume accordingly: Assumption 2.1 The matrix function P = B(GB) 1G is smooth, or more precisely, P is continuous and P 0 is bounded near the singular point t where we de ne P (t ) = lim t!t (B(GB) 1G)(t): Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations23 Because we are only interested in the case where (2.16) has a smooth solution for x (as is the case in Example 2.1), it is necessary to assume, in view of (2.16b): Assumption 2.2 The inhomogeneity r(t) satis es r 2 S; where S = fw(t) 2 Rny : Gz = w for a smooth function z(t): We note that Assumptions 2.1 and 2.2 are satis ed automatically if GB is nonsingular for each t. On the other hand, neither B(GB) 1 nor (GB) 1G alone are smooth near a singularity in general. We also indicate here that to formulate the SRM (see x2.4) we only need the continuity of P . The further requirement in Assumption 2.1 on the derivative of P is needed for the regularity of the solution (cf. Lemma 2.1 and (2.24)) and the stability proof for the regularized problems (cf. x2.5). This requirement can be avoided if we make a more general assumption about the regularity and stability of the original problem (cf. x3.1). We consider both initial and boundary value problems. In x2.3 we brie y discuss the conditioning of the problem (2.16) with singularities. In x2.4 we derive the sequential regularization method. In x2.5 we estimate the error of the SRM. In x2.6, we consider some discretization and implementation issues for both initial and boundary value problems. Finally, in x2.7 several numerical examples demonstrate our theoretical results. 2.3 Problem Conditioning Similarly to [15] and to the method of pseudo upper triangular decomposition (PUTD) described in [2] (cf. x10.6; with the di erence that we do pivoting to interchange the row with the singularity of the lowest order and the current row when all the Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations24 other rows vanish at some singular point), there exists a smooth matrix function R(t) 2 R(nx ny) nx , which has full row rank and satis es RB = 0; for each t; 0 t tf ; where R can be taken to have orthonormal rows. As in [15, 16], de ne the new variable u = Rx; 0 t tf : (2.19) Then, using (2.16b), the inverse transformation is given by x = Su B(GB) 1r; (2.20) where S = (I B(GB) 1G)RT = (I P )RT : By the assumptions at the beginning of this chapter, this transformation is wellde ned. Di erentiating (2.19) and using (2.16a) and (2.20) we obtain the essential underlying ODE (EUODE): u0 = (RA+R0)Su (RA+R0)B(GB) 1r +Rq: (2.21) Hence the underlying problem of (2.21) is u0 = (RA+R0)Su+ f; (2.22a) B0S(0)u(0) + B1S(tf)u(tf) = 1: (2.22b) We make Assumption 2.3 The boundary value problem (2.22) is stable, i.e. there exists a moderate-size constant K such that kuk K(kfk+ j 1j); where kuk = maxtfju(t)j; 0 t tfg. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations25 Similarly to Theorem 2.2 of [15], we have Lemma 2.1 Let the DAE (2.16) have smooth coe cients, and assume that Assumptions 2.1 and 2.2 hold. If the EUODE (2.21) with boundary conditions (2.22b) has a unique solution, then there exists a unique solution for the x of problem (2.16)(2.17) which is smooth. This implies the unique existence of a smooth By as well. Furthermore, if Assumption 2.3 holds then there is a constant K such that kxk K(kqk+ kB(GB) 1rk+ j j); kx0k K(kqk+ kB(GB) 1rk+ k(B(GB) 1r)0k+ j j): Remark 2.1 For problem (2.16) without singularities, we can get kxk K(kqk+ krk+ j j); kx0k K(kqk+ krk+ kr0k+ j j) as in [15, 16]. 2 The di erence between the situation here and in the nonsingular case is that the perturbation inhomogeneities r yield reasonably bounded perturbations in the solution x only if they are (in general) from the subspace Range (G). From (2.16a) and (2.20), we can write y = (GB) 1G(x0 Ax q); t 2 [0; t ) [ (t ; tf ]; (2.23) which could be unbounded at the singular point t (whereas By is bounded). Note that G could be replaced in (2.23) by any appropriate matrix Q with the same size as G, e.g. Q can be BT . Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations26 Remark 2.2 If B has full rank for each t, then (GB) 1G = (BTB) 1BTP: Hence, (GB) 1G is smooth. Hence, there exists a unique solution for the y of problem (2.16)-(2.17) which is smooth and can be expressed as (2.23) for each t. Furthermore, using Lemma 2.1, we have in this case kyk K(kqk+ kB(GB) 1rk + k(B(GB) 1r)0k+ j j): (2.24) In the general case, however, we will have to consider By, rather than y alone, in the theorems of the next section. 2 A Baumgarte stabilization applied to (2.16) consists of eliminating y according to (2.18),(2.23), and stabilizing (see (2.30) below for the usual form). This gives the ODE x0 = (I B(GB) 1G)(Ax+ q) B(GB) 1(G0x+ r0) 1B(GB) 1(Gx+ r) (2.25) where > 0 is a parameter (cf. [17, 8]). If there are no singularities then it follows from the analysis in [16] that if Assumption 2.3 holds then the boundary value problem (2.25),(2.17),(2.18b) is also stable. In other words, the \initial value stabilization" works also for the boundary value case, because the new modes introduced by replacing (2.16b) with (2.18a) are separable and decaying, in agreement with the additional initial conditions (2.18b). However, in the singular case (2.25) may not work because the terms B(GB) 1G0 and B(GB) 1r0 are in general unbounded. Therefore, we develop an iterative method in the next section which builds up an approximation to By and x that avoids going through unbounded quantities. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations27 2.4 Derivation of the SRM There are many discussions on regularization methods for DAEs. A direct regularization (cf. the pseudo-compressibilitymethod in the Navier-Stokes context [108, 72, 97]) is: x0 = A(t)x B(t)y + q(t); (2.26a) y0 = G(t)x+ r(t): (2.26b) This formulation is not popular because it requires conditions on A, B and G, for the purpose of stability of the system, and the existence of the rst derivative of y, which is not necessarily true for the original problem (2.16) [86]. It may also change the properties of the original index-2 problem too much by jumping from index-2 to index-0 (ODE). It seems evident that a regularization with fewer changes of the original problem (e.g. from index-2 to index-1) might be better. The penalty method [84, 91, 68, 70] is such a method. It reads x0 = A(t)x B(t)y + q(t); (2.27a) E 1y = G(t)x+ r(t); (2.27b) where E 2 Rny ny is chosen such that BEG has non-negative eigenvalues. Hence, the system obtained by substituting (2.27b) into (2.27a) is generally stable. For example, we can choose, relying on Assumption 2.1, E = (GB) 1(hence, BEG = P ). Also, E = (GB)T could be a good choice in some circumstances. If B =M 1GT for some positive de nite matrix M (as in the case of mechanical systems) then it is possible to choose E = I. Advantages of these choices of E will be discussed in Chapter 3. For problems with singularities, we suggest using E = (GB) 1 to avoid a turning point problem. Another approach is parameterization [85, 59]: x0 = A(t)x B(t)y+ q(t); (2.28a) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations28 0 = G(t)(x+ x0) + r(t): (2.28b) Substituting (2.28a) to (2.28b), we get GBy = Gx+ r + GAx: (2.29) This implies that parameterization can be seen as an instance of the penalty method with E = (GB) 1. Recently, [94] has reported a regularization for DAEs with singularities based on a formulation obtained by the trust-region method in numerical optimization. All these regularizations require the regularization parameter to be very small. Therefore a sti solver is needed to solve the regularized problem. In this section, we derive a new regularization method which is called the sequential regularization method [11]. The SRM is an iterative procedure which combines the popular Baumgarte stabilization or other stabilizations with a modi ed penalty method. One purpose for doing so is that the regularized problems of the SRM can be non-sti or at least less sti . Hence, simple discrete schemes (e.g. explicit schemes) can be used. Other advantages of the method will be discussed in Chapter 3. The Baumgarte stabilization of (2.16) reads (cf. (2.5)) 1 d dtg(x; t) + 2g(x; t) = 0; g(x(0); 0) = 0: (2.30) Applying the idea of the penalty method to equations (2.16a) with constraints (2.30), we obtain x0 = A(t)x B(t)y + q(t); (2.31a) y = y0 + 1 E( 1 d dtg(x; t) + 2g(x; t)): (2.31b) where y0 can be seen as an initial guess of the exact solution ye of problem (2.16), (2.17). If we take y0 = ye then the solution of problem (2.31), (2.17) is exactly the Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations29 same as that of problem (2.16), (2.17). If we take y0 0 then (2.31) coincides with the penalty method (2.27). Given any initial guess y0(t), the solution, say fx1; y1g, of (2.31), (2.17) is an approximation of the exact solution fxe; yeg of (2.16), (2.17). Using this solution y1 as a new initial guess, we re-solve problem (2.31), (2.17). We expect that the solution obtained is a better approximation of the exact solution. Repeating the procedure, we invent the following iterative algorithm for solving (2.16): For s = 1; 2; : : : , solve the ODE problem xs0 = Axs Bys + q (2.32) where ys = ys 1 + 1 E( 1 d dtg(xs; t) + 2g(xs; t)); (2.33) subject to the same boundary conditions (2.17). Note that y0(t) is a given initial iterate and that > 0 is the regularization parameter. We call this algorithm a sequential regularization method (SRM). Note that xs(t) and ys(t) are de ned on the entire interval [0; tf ] for each s. For the problem with singularities, the choice E = (GB)T generates turning point regularized problems which are complicated to solve and analyze. We thus choose E = (GB) 1 for the singular case. Noting that B(GB) 1G0 may be unbounded at the singularity we then choose 1 = 0 to avoid this term. Also, in practice we multiply (2.33) by B and keep track only of the approximations Bys to the bounded function By, since y may be unbounded at the singularity. We thus have an SRM variant for the singular case: For s = 1; 2; : : : , solve the ODE problem xs0 = Axs Bys + q (2.34) where Bys = Bys 1 + 1 BEg(xs; t); (2.35) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations30 subject to the boundary conditions (2.17). If y is desired (at times other than at the singular point t ) then it can be easily retrieved from By in a post-processing step. 2.5 Convergence Analysis of the SRM We rst prove a lemma which will be used to discuss the convergence of the SRM. Lemma 2.2 Let u, v be the solution of u0 = (RA+R0)Su+ S1v + f1; (2.36a) v0 + v = S2u+ S3v + f2; (2.36b) B0S(0)u(0) + B1S(tf)u(tf) = S4v(0) S5v(tf); v(0) = v0; (2.36c) where all coe cients are bounded, = 1 or = , is a positive constant and Assumption 2.3 holds. Then, for appropriately small or appropriately large, we have the following stability inequality kuk K(kf1k+ kf2k+ j j+ jv0j); kvk K( kf1k+ kf2k+ j j+ jv0j); where K is a positive constant. Proof: Let v = (v1; ; vny)T . From (2.36b), we easily have jvij kS2kkuk+ kS3kkvk+ 1 kf2k+ jv0j; i = 1; ; ny: Hence, taking the maximum of the left hand side for 1 i ny and choosing small or large appropriately such that kS3k < , we get kvk kS3kkS2kkuk+ kf2k+ jv0j kS3k : (2.37) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations31 By using Assumption 2.3, from (2.36a), there exists a positive constant K1 such that kuk K1(kS1kkvk+ kf1k+ j j+ jS4jjv(0)j+ jS5jjv(tfj) K1((kS1k+ jS5j)kvk+ kf1k+ j j+ jS4jjv0j) K1 (kS1k+ jS5j)kS2k kS3k kuk+K1(kS1k+ jS5j)(kf2k+ jv0j) kS3k +K1(kf1k+j j+jS4jjv0j): Hence, by choosing smaller or larger such that K1 (kS1k+jS5j)kS2k kS3k < 1, the stability inequality for u follows. Now the stability inequality for v follows from that for u and (2.37). 2 Now we estimate the error of the SRM (2.34), (2.35). De nition 2.1 J is an integer such that y0(0) = ye(0); y0 0(0) = y0 e(0); : : : ; y(J) 0 (0) = y(J) e (0); where y0(t) is the initial guess of the SRM iteration (2.34), (2.35) and ye(t) is the exact solution of the original problem (2.16). Set J = 1 if y0(0) 6= ye(0). 2 For initial value problems we may calculate ye(i)(0); i = 0; 1; : : : in advance by using the ODE and its derivatives. For boundary value problems we have J = 1 in general since we don't know ye(0) beforehand. Theorem 2.1 Let the DAE (2.16) have su ciently smooth coe cients, and assume that Assumptions 2.1, 2.2 and 2.3 hold. Then, for the solution of iteration (2.34),(2.35), we have the following error estimates ( for J de ned in De nition 2.1): xs(t) xe(t) = O( s) +O( J+2ps(t= )e t= ); Bys(t) Bye(t) = O( s) +O( J+1ps(t= )e t= ); for 0 t tf and s 1. Here ps( ) 0 if s J+1; otherwise ps( ) is a polynomial of degree s J 2 with generic positive coe cients and jps(0)j = j(By0)(J+1)(0) (Bye)(J+1)(0)j. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations32 Proof: Let us = Rxs and ws = Pxs. Similarly to (2.20), we have xs = Sus + ws: (2.38) Furthermore, using (2.34) we obtain u0s = (RA +R0)Sus + (RA+R0)ws +Rq; (2.39a) w0 s + ws = (PA+ P 0)Sus + (PA+ P 0)ws Bys 1 (2.39b) + Pq B(GB) 1r; subject to B0S(0)us(0) + B1S(tf)us(tf ) = B0ws(0) B1ws(tf); (2.40a) ws(0) = B(0)(G(0)B(0)) 1r(0): (2.40b) The iteration (2.35) for By becomes Bys = Bys 1 + 1 (ws +B(GB) 1r): (2.41) The proof proceeds along familiar lines of singular perturbation analysis. According to [111, 112] we can construct the asymptotic expansion of ws and us sequentially for s = 1; 2; : : :, where we use Lemma 2.2 to estimate the remainders. Then, using (2.41) and (2.38), we get the asymptotic expansion of Bys and xs respectively. Note that in these expansions the rst terms are exactly xe and Bye. This process eventually yields the proof of the theorem. 2 To provide a better understanding about the sequential regularization method we give in x2.8 a detailed proof for the initial value case with no layers, s J + 1. In that proof, the construction of the asymptotic expansion is directly for x and By. Moreover, the construction method we apply is somewhat di erent from [111, 112] and more relevant to the concept of DAEs. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations33 Next, we consider the SRM (2.32), (2.33). For the initial value problem with E = I and B = GT , this corresponds to Algorithm (2.13), (2.14) of [19] for constrained mechanical systems (although they do it for the corresponding index-3 case) derived by a penalty-augmented Lagrangian formulation. Bayo and Avello have noted that under repetitive singular conditions this algorithm may lead to unstable behavior. For our index-2 case (2.16) with singularity, it appears to be impossible to choose a matrix E such that problem (2.32),(2.33) is always stable, even if we assume B = GT . A numerical example in x2.7 will verify such instability phenomena even for the case of one singular point. However, for the case where constraints are without singularities, (2.33) (multiplied by B) may have a bene t over (2.35). That is, (2.33) yields an ODE problem for xs which is essentially not a sti problem. Take E = (GB) 1 as before and rewrite (2.33) as Bys = Bys 1 + 1 BE( 1 d dtg(xs; t) + 2g(xs; t)): (2.42) Then we give the following error estimation for (2.32), (2.42): Theorem 2.2 Let the DAE (2.16) have su ciently smooth coe cients, and assume that G has full rank and that Assumptions 2.1, 2.2 and 2.3 hold. Then for the solution of the iteration procedure (2.32), (2.42) with 1 6= 0, we have the following error estimates: xs xe = O( s); Bys Bye = O( s) for 0 t tf and s = 1; 2; : : :. Note that no boundary layer terms appear here even for J = 1 (See De nition 2.1)! Proof: Denote us = Rxs and vs = Gxs. Hence xs = Sus + Fvs; (2.43) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations34 where S = (I P )RT and F = B(GB) 1 = PGT (GGT ) 1 are both su ciently smooth. From (2.32),(2.42), we get u0s = (RA+R0)Sus + (RA+R0)Fvs +Rq; ( + 1)v0 s + 2vs = (G0 +GA)Sus + (G0 +GA)Fvs + GBys 1 + Gq 1r0 2r; with the corresponding boundary conditions, and Bys = Bys 1 + 1 B(GB) 1( 1(vs + r)0 + 2(vs + r)): Repeating the procedure of the proof of Theorem 2.1 and using Lemma 2.2 again to estimate the remainder of the asymptotic expansion, we obtain us ue = O( s); vs ve = O( s); Bys Bye = O( s); where ue = Rxe; ve = Gxe = r. Hence, using (2.43) and xe = Sue+Fve, we obtain xs xe = S(us ue) + F (vs ve) = O( s): 2Remark 2.3 For the problem (2.32), (2.33) where GB is nonsingular, we have xs xe = O( s); ys ye = O( s): This estimate also holds for E = (GB)T . 2 2.6 Discretization and Implementation Issues The SRM iteration yields a sequence of ODE problems which are to be solved numerically. We only consider the most di cult case, i.e. (2.34), (2.35) with singularities. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations35 Inserting (2.35) into (2.34), the ODE problem to be solved at the sth iteration is written as the singular-singularly-perturbed problem (see [112, 89]) x0s +BE(Gxs + r) = Axs (Bys 1 q); (2.44a) B0xs(0) + B1xs(tf ) = ; G(0)xs(0) + r(0) = 0: (2.44b) We consider nite di erence (or collocation) discretizations of (2.44) on a mesh : 0 = t0 < t1 < : : : < tN = tf hi = ti ti 1; h = max 1 i N hi and denote by xsi , ys i the corresponding approximations of xs(ti), ys(ti), respectively. We now have essentially two small, positive parameters to choose: and h. We assume that h is chosen small enough so that the EUODE problem (2.22) may be considered as nonsti and that the problem's coe cients are su ciently smooth. In the BVP case the situation is the familiar one, much like the iterative solution of a nonlinear boundary value ODE using quasilinearization (see, e.g., [14]). Each of the linear boundary value ODEs (2.44) is discretized on a mesh using, say, a symmetric nite di erence scheme or some other method. We expect, as h! 0, convergence to the solution of ( 2.44) and our theory then applies for the entire numerical algorithm. As an example, we give here a detailed analysis of the convergence of the backward Euler di erence scheme for (2.44). A similar discussion and results can easily apply to the forward Euler di erence scheme. The results for general higher order collocation schemes have been described in [11]. Now we write the backward Euler scheme of (2.44) as follows: xsi xsi 1 hi +BiEi(Gixsi + ri) = Aixsi (Biys 1(ti) qi); (2.45a) B0xs0 + B1xsN = ; G0xs0 + r0 = 0; (2.45b) Biys i = Biys 1(ti) + 1 BiEi(Gixsi + ri); (2.46) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations36 where we represent the value f(ti) of a given function f at mesh point ti by fi. Multiplying (2.45) by Ri and Pi, respectively, and denoting usi = Rixsi ; ws i = Pixsi (then xsi = Siusi + ws i since (2.38) holds), we have usi usi 1 hi = RiAi(Siusi + ws i ) + Ri Ri 1 hi (Si 1usi 1 + ws i 1) +Riqi; (2.47a) ws i ws i 1 hi + ws i = PiAi(Siusi + ws i ) + Pi Pi 1 hi (Si 1usi 1 + ws i 1) (Biys 1(ti) qi) BiEiri; i = 1; ; N; (2.47b) B0(S0us0 + ws 0) + B1(SNusN + ws N) = ; ws 0 = B0E0r0 = P0xe(t0): (2.47c) At rst, we consider the stability of the following di erence scheme corresponding to (2.47): Lh1ui ui ui 1 hi RiAiSiui Ri Ri 1 hi Si 1ui 1 = RiAiwi + Ri Ri 1 hi wi 1 + fi; (2.48a) Lh2wi wi wi 1 hi + wi = PiAiSiui + Pi Pi 1 hi Si 1ui 1 PiAiwi + Pi Pi 1 hi wi 1 + gi; i = 1; ; N; (2.48b) B0S0u0 + B1SNuN = 1 = B0w0 B1wN ; w0 = 2: (2.48c) Using the discrete maximum principle for the di erence operator Lh2, i.e. z0 0 and Lh2zi 0 =) zi 0;8i; (2.49) we easily get max 1 j i jwjj M(j 2j+ max 1 j i jLh2wjj): (2.50) De ning kzki = max1 j i jzjj and using (2.48b), we have kwki M (kuki + kwki) +M j 2j+ kgki Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations37 or kwki M( kuki + j 2j+ kgki): (2.51) Here M stands for a generic constant independent of i, and h. On the other hand, Lh1ui is a one-step di erence operator of (2.22a) | the underlying problem of (2.21). Using Assumption 2.3 and Theorem 5.38 of [14], we obtain kuk1 K1(j 1j+ max 1 j N jLh1ujj); (2.52) where kzk1 = max0 j N jzjj and K1 = K + O(h). Here K is the stability constant de ned in Assumption 2.3. Using (2.48a) and (2.48c) we have kuk1 M(kwkN + j 2j+ j j+ kfkN ): (2.53) Then, using (2.51), yields kuk1 M(kfkN + kgkN + k j+ j 2j) (2.54a) kwk1 M( kfkN + kgkN + k j+ j 2j): (2.54b) We thus obtain the stability inequalities (2.54) for the di erence scheme (2.48) or (2.47). Now we discuss the convergence of (2.47). Using (2.54) , we have the estimates: kus(ti) usik1 M(k ukN + k wkN ) (2.55a) kws(ti) ws i k1 M(k ukN + k wkN ); (2.55b) where u i and w i are local truncation errors for the di erence scheme (2.47) and they can be written as u i = hifu00 s( 1 i ) +R00( 2 i )[S(ti 1)us(ti 1) + ws(ti 1)] + R0(ti)(Sus + ws)0( 3 i )g (2.56a) w i = hifw00 s ( 1 i ) + P 00( 2 i )[S(ti 1)us(ti 1) + ws(ti 1)] + P 0(ti)(Sus + ws)0( 3 i )g; (2.56b) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations38 where i ; i 2 (ti 1; ti); = 1; 2; 3. To bound the truncation error, we need the derivative estimates of us and ws. From the asymptotic expansions of us and ws ( cf. Theorem 2.1) us = ue + us1 + 2( us2 + ~ us2) + (2.57a) ws = we + ( ws1 + ~ ws1) + 2( ws2 + ~ ws2) + (2.57b) where ue = Rxe; we = Pxe , usj and wsj are functions of regular expansions, ~ usj and ~ wsj are boundary layer functions whose basic forms are p(t= ) exp( t= ) (where p is some polynomial) and usj = wsj = 0 for j s ( since SRM iteration cancels out the lower terms of regular expansions). We can expect that ju0sj; ju00 sj; jw0 sj M: (2.58) But jw00 s j = O( 1 exp( t= ) 8<: M if t M 1 otherwise : (2.59) Therefore k ukN = O(h); k wk = 8<: O( h) if h1 = t1 t0 O(h) otherwise : (2.60) From (2.55), we thus have kus(ti) usik1 Mh (2.61a) kws(ti) ws i k1 8<: M h if h1 Mh otherwise (2.61b) i.e. xs(ti) xsi = S(ti)(us(ti) usi ) + (ws(ti) ws i ) = O(h): (2.62) However, we can not generally get a good approximation for Bye in the whole region if is not very small compared with h1 since in this case we generally have Biys i = Biys 1(ti) + 1 BiEi(Gixsi + ri) = Biys 1(ti) + 1 BiEi(Gixs(ti) + ri) + 1 (ws i ws(ti)) (2.63) = Biye(ti) +O( + exp( ti= )) +O(h= ): Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations39 Fortunately, we can get O(h) accuracy locally, i.e. in a smooth region or away from the layer region, say for ti0 t tf . Indeed, considering an equidistant mesh for simplicity, from (2.47b), we have Lh1 ws i = ws i ws i 1 h + ws i 1 = PiAi(Si usi + ws i ) + Pi Pi 1 h (Si 1 usi + ws i 1) +O( w i ); (2.64) where ws i = ws(ti) ws i , usi = us(ti) usi and we note that w i = O( h) for i0 i N . Using the discrete maximum principle for Lh1 in ti0 t tf on the barrier function zi = j ws i0j i i0 + max i0 i N j ij ws i ; where i(= O(h)) is the right-hand side of (2.64) and = h=( +h), and using (2.61), we get zi 0 or j ws i j j ws i0j i i0 + max i0 i N j ij M(h i i0 + h): (2.65) For = h1+ ; 1 < < 1, i i0 = exp( (i i0)h ) = exp( (ti ti0)=h1 ): Taking i1 such that exp( (ti1 ti0)=h1 ) Mh1+ = M when h is su ciently small, we get i i0 exp( (ti1 ti0)=h1 ) M for i i1; i.e. j ws i j M h; 8i i1, and any satisfying = h1+ ; 1 < < 1. Combining this with (2.61b), we obtain j ws i j M h; 8i i1: (2.66) Then, using (2.63) and (2.66), we have the local error estimate jBiys i Biye(ti)j M(h + s); for i1 i N: (2.67) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations40 This means we can get good approximation in a region which is away from the initial layer. Once an accurate SRM solution, say fx i ; B(ti)y i g, has been determined outside the initial layer it may be possible to obtain an accurate solution everywhere by applying a few SRM iterations numerically solving (2.44a) ( changing BEG to BEG) subject to the terminal value x(tf) = x N ; (2.68) and choosing By0 satisfying B(tf)y0(tf) = B(tf)y N . This procedure is feasible provided that the terminal value problem (2.44a),(2.68) is well-conditioned (which holds if the terminal value problem for the EUODE (2.22a) is well-conditioned). For the IVP case, where (2.44b) reduces to x(0) = x given; (2.69) we may, of course, proceed in the same way as for the BVP case. But now a few things are easier. Firstly, for this case we can calculate Bye(0) and then choose By0 to be exact at t = 0. In fact, as indicated earlier we can also do this for higher derivatives of By at the initial value by repeated di erentiation of (2.16). Such a preparation of the initial iterate By0 allows removing the layer error terms (or the condition ti ) in the error estimates (2.61). Secondly, one may use a more convenient di erence scheme to integrate the IVP (2.44a),(2.69). If the EUODE is su ciently nonsti to warrant use of a nonsti integration method then this can be an attractive possibility here. Note, though, that hi= must be in the absolute stability region of the method (see (2.39b)). Thus, an explicit Runge-Kutta method of order p, for instance, may necessitate (at least) p SRM iterations in order for the error in the estimates of Theorem 2.1 to be of the same order as the error in the numerical approximation. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations41 The most important di erence between the IVP and BVP cases is that the iterative method described here does not appear to be necessarily optimal or even natural in the IVP context, certainly not from the storage requirement point of view: Note that the entire approximation of Bys 1 on [0; tf ] needs to be stored. The situation here is similar to that encountered with other functional iteration methods like waveform methods. However, this di culty can be resolved by rearranging the computation, assuming that the number of the SRM iterations, m, is chosen in advance. Thus, at each time step i, 1 i N , we calculate sequentially the quantities x1i ; By1 i ; x2i ; By2 i ; : : : ; xmi ; Bym i . To do this using a one-step scheme, say, we need only the corresponding quantities locally, over the mesh subinterval [ti 1; ti), and By0 i . For the latter we may use, for instance, y0 i ye(0), i.e. By0 i = B(ti)y0 0; 0 i N . The storage requirements are now independent of N and other typical IVP techniques such as local error control may be applied as well. 2.7 Numerical Experiments We now present a few very simple examples to demonstrate our claims in the previous sections. Throughout this section we use a constant step size h and set tf = 1. To make life di cult we choose h so that there is an i such that ti = t (if there is a singularity). In the implementation we monitor the size of the pivot in a Gauss elimination procedure for GB and slightly perturb ti away from t when needed. At a given time t, we use 0ex0 to denote the maximum over all components of the error in xs while 0ey0 denotes the maximum over all components of the error in Bys. Similarly, 0drift0 denotes the maximum residual in the algebraic equations. We rst look at a boundary value problem. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations42 Example 2.2 Consider the DAE (2.16) with A = 0@ 1 1 0 01A ; B = 0@ 0 1 2t1A ; q = 0@ sin t 0 1A G = ( 1 2t 1 2t ) ; r = (1 2t)(e t + sin t) subject to x1(1) + x2(0) = 1=e: The exact solution is xe = ( e t sin t ), ye = cos t 1 2t . A singularity is located at t = 1=2, where ye becomes in nite while Bye stays bounded. We start computing with the iterate y0(t) 0. In Table 2.1 we list errors when using the midpoint scheme xsi xsi 1 hi = Ai 1 2xsi 1 2 Bys i 1 2 + qi 12 Bys i 12 = Bys 1 i 12 + 1Bi 12Ei 12 (Gi 12xsi 12 + ri 1 2 ) where xsi 1 2 = xsi+xsi 1 2 (but no such relation is necessary for ys). We apply this scheme with hi = h = :01 for various values of . It is indicated in [11] that this scheme has 2nd order accuracy in ex and in ey, except for the case h when the error's order in By drops to 1. This is evident in the error column for t = 1:0. Note also the O( ) improvement per SRM iteration when this term dominates the error (i.e. when s h2). We note that the approximation for By at points within the initial layer is not accurate. To get a better approximation within the initial layer ( i.e. near the initial point t = 0), we solve a terminal value problem (2.44a) (changing BEG to BEG), (2.68), as described in x2.6. Then we apply the SRM for the given problem with the improved values for By0. In Table 2.2 we list the computed results after 3 iterations. They are obviously much better than the comparable ones in Table 2.1. 2 Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations43 iteration error at ! t=.01 t=.1 t = .3 t=.5 t=1.0 1e-1 1 ex .38e-1 .35e-1 .56e-1 .52e-1 .39e-1 ey .96 .40 .14 .34e-1 .66e-1 drift .87e-2 .49e-1 .51e-1 .0 .61e-1 2 ex .92e-2 .37e-1 .89e-2 .65e-2 .72e-2 ey .91 .96e-2 .14 .34e-1 .65e-2 drift .90e-2 .32e-1 .61e-2 .0 .72e-2 3 ex .94e-2 .19e-1 .12e-1 .63e-2 .15e-2 ey .87 .20 .30e-1 .43e-1 .85e-3 drift .86e-2 .16e-1 .43e-2 .0 .74e-3 1e-2 1 ex .38e-2 .60e-2 .53e-2 .44e-2 .38e-2 ey .67 .30e-2 .40e-2 .48e-2 .62e-2 drift .65e-2 .80e-2 .38e-2 .0 .55e-2 2 ex .45e-2 .10e-3 .88e-4 .77e-4 .64e-4 ey .45 .32e-3 .70e-4 .44e-4 .23e-4 drift .44e-2 .15e-4 .14e-4 .0 .68e-4 3 ex .30e-2 .55e-5 .52e-5 .59e-5 .11e-4 ey .30 .18e-2 .26e-4 .18e-4 .67e-5 drift .29e-2 .17e-5 .26e-5 .0 .56e-5 1e-3 1 ex .13e-2 .58e-3 .52e-3 .45e-3 .39e-3 ey .17 .47e-2 .41e-3 .49e-3 .62e-3 drift .17e-2 .79e-3 .38e-3 .0 .54e-3 2 ex .30e-3 .71e-4 .75e-5 .72e-5 .12e-4 ey .30e-1 .17e-1 .34e-4 .13e-4 .54e-5 drift .30e-3 .51e-4 .20e-5 .0 .65e-5 3 ex .65e-4 .15e-3 .70e-5 .70e-5 .12e-4 ey .70e-2 .33e-1 .12e-3 .14e-4 .56e-5 drift .69e-4 .11e-3 .21e-5 .0 .59e-5 Table 2.1: SRM errors for Example 2.2 using the midpoint scheme Next we consider initial value problems. Example 2.3 Consider the same DAE as for Example 2.2 with the same exact solution but with initial values x1(0) = 1, x2(0) = 0 speci ed. From these initial conditions we can calculate y(0) = 1 in advance, and we choose the initial guess y0(t) 1. Tables 2.3 and 2.4 display error results for = :1 and h = :001 using the backward Euler and the forward Euler schemes, respectively. As explained in x2.6 we calculate all iterates at each step before proceeding to the next. These tables show a signi cant improvement with each SRM iteration and no strong initial layer e ect, as predicted by theory. 2 Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations44 error at ! t=.01 t=.1 t = .3 t=.5 1e-1 ex .62e-2 .54e-2 .39e-2 .28e-2 ey .48e-2 .45e-2 .44e-2 .54e-2 5e-2 ex .76e-3 .58e-3 .32e-3 .17e-3 ey .22e-2 .18e-2 .10e-2 .50e-3 1e-2 ex .57e-4 .49e-4 .35e-4 .26e-4 ey .10e-3 .85e-4 .52e-4 .32e-4 1e-3 ex .49e-4 .42e-4 .31e-4 .23e-4 ey .56e-4 .46e-4 .30e-4 .19e-4 Table 2.2: SRM errors for Example 2.2 using the shooting-back technique iteration error at ! t=.001 t=.1 t = .3 t=.5 t=1.0 1 ex .20e-5 .72e-2 .37e-1 .63e-1 .11 ey .20e-2 .12 .15 .12 .59e-1 drift .15e-5 .60e-2 .16e-1 .76e-4 .15 2 ex .20e-5 .51e-2 .13e-1 .10e-1 .25e-2 ey .20e-2 .68e-1 .45e-2 .20e-1 .80e-2 drift .15e-5 .42e-2 .58e-2 .14e-4 .67e-2 3 ex .20e-5 .35e-2 .23e-2 .16e-2 .76e-3 ey .20e-2 .32e-1 .26e-1 .10e-1 .37e-2 drift .15e-5 .29e-2 .12e-2 .10e-5 .12e-2 Table 2.3: SRM errors for Example 2.3 using backward Euler Example 2.4 Here we investigate the use of the modi ed formula (2.42) instead of (2.35). First, we solve the previous example numerically using (2.42). In Table 2.5 we record error values at the singularity point t = :5 after 3 SRM iterations, starting with y0(t) 1 and using as before = :1 and h = :001 (cf. Tables 2.5). From these results it is clear that the SRM with (2.42) does not work well when 1 6= 0: large errors in By are obtained near the singularity and these adversely a ect the accuracy in x as well. However, the comparison changes when there is no singularity in the constraints: We now replace the algebraic constraint in Example 2.3 by x1 + x2 e t sin t = 0 leaving everything else the same (including the singularity in B). In Table 2.6 we record maximum errors in x and By over all mesh points (denote those 'exg' and Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations45 iteration error at ! t=.001 t=.1 t = .3 t=.5 t=1.0 1 ex .50e-6 .71e-2 .36e-1 .63e-1 .11 ey .20e-2 .12 .15 .12 .60e-1 drift .50e-6 .60e-2 .16e-1 .76e-4 .15 2 ex .50e-6 .51e-2 .12e-1 .10e-1 .44e-2 ey .20e-2 .68e-1 .41e-2 .20e-1 .70e-2 drift .50e-6 .42e-2 .58e-2 .14e-4 .67e-2 3 ex .50e-6 .35e-2 .43e-2 .18e-2 .98e-3 ey .20e-2 .32e-1 .26e-1 .97e-2 .46e-2 drift .50e-6 .29e-2 .12e-2 .11e-5 .12e-2 Table 2.4: SRM errors for Example 2.3 using forward Euler ( 1; 2) ! (0; 1) (h; 1) (1; 1) method ex ey ex ey ex ey backward Euler .16e-2 .10e-1 .80e-3 .20e+1 .15 .37e+3 forward Euler .18e-2 .96e-2 .25e-2 .18e+1 .64 .15e+4 Table 2.5: Errors near singularity using modi ed formula (2.42) 'eyg', respectively) for the starting iterates y0(t) 1 and y0(t) 0 (the latter does not agree with the exact ye(0)). ( 1; 2) ! (0; 1) (h; 1) (1; 1) y0 method exg eyg exg eyg exg eyg 1 backward Euler .46e-2 .44e-1 .45e-2 .43e-1 .22e-3 .28e-3 forward Euler .45e-2 .44e-1 .44e-2 .43e-1 .19e-3 .23e-3 0 backward Euler .22e-1 .97 .22e-1 .94 .12e-3 .75e-3 forward Euler .22e-1 .97 .22e-1 .94 .19e-3 .75e-3 Table 2.6: Errors for problem without singularity using modi ed formula (2.42) The modi ed method (corresponding to ( 1; 2) = (1; 1) in Table 2.6 is seen to work better for problems without singularities. 2 The above calculations all agree with our theoretical results described in x2.5 and x2.6. Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations46 2.8 More about the Proof of Theorem 2.1 To provide a better understanding about the sequential regularization method we now give a detailed proof of Theorem 2.1 for the initial value case with no layers, s J+1. J is some positive integer de ned in De nition 2.1. In this proof, the construction of the asymptotic expansion is directly for x and By. Moreover, the construction method we apply is somewhat di erent from [111, 112] and more relevant to the concept of DAEs. The same idea is applied to prove the convergence of the SRM for NavierStokes equations in Chapter 4. For s > J+1, additional initial layer expansions have to be developed. However, the construction of these layer expansions is precisely the same as in [111, 112] and so it is omitted here. In case that (2.17) are initial conditions (i.e. B1 = 0) our assumptions imply that (2.17) together with (2.18b) specify x(0), say x(0) = x (2.70) At rst, consider the case s = 1 of (2.34),(2.35): x01 +B(GB) 1(Gx1 + r) = Ax1 By0 + q; with the initial conditions (2.70). This is a singular-singularly-perturbed problem (see [112, 89]). Let x1 = x10 + x11+ + sx1s + Comparing the coe cients of like powers of , we thus have B(GB) 1Gx10 = B(GB) 1r (2.71a) B(GB) 1Gx11 = x010 +Ax10 By0 + q; (2.71b) B(GB) 1Gx1i = x01i 1 +Ax1i 1; 2 i s+ 1; (2.71c) where (2.71a) satis es (2.70) and (2.71b) and (2.71c) satisfy homogeneous initial conditions corresponding to (2.70). Now, (2.71a) has in nitely many solutions in Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations47 general. To realize the construction, we should choose x10 to satisfy (2.71a) and to ensure that the solution of (2.71b) exists. We choose x10 to be the solution xe of problem (2.16)-(2.17), i.e. x010 = Ax10 Bye + q; (2.72a) 0 = Gx10 + r; (2.72b) B0x10(0) = : (2.72c) So x10 = xe and (2.71b) has the following form B(GB) 1Gx11 = B(ye y0): (2.73) Now we choose x11 and a corresponding y01 to satisfy x011 = Ax11 By01 (2.74a) Gx11 = GB(ye y0); (2.74b) B0x11(0) = 0: (2.74c) Noting that Bye = x0e + Axe + q is smooth, we have GB(ye y0) 2 S. Hence, using Lemma 2.1, there exists a smooth solution x11 of (2.74), and x11 satis es (2.73). Indeed, using (2.74b) and De nition 2.1, we have G(0)x11(0) = 0, so x11(0) = 0. And, from (2.74b) again, (GB) 1Gx11 = ye y0; for each t 2 [0; t ) [ (t ; tf ]: That is, B(GB) 1Gx11 = B(ye y0); t 2 [0; t ) [ (t ; tf ]: (2.75) Taking the limit of (2.75) at t , we thus get that x11 satis es (2.73) for each t 2 [0; tf ]. Moreover, from De nition 2.1, we have y01(0) = y0 01(0) = = y(s 1) 01 (0) = 0; s J + 1: Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations48 Also we note that By01 is smooth. Generally, supposing we have got x1i 1; By0i 1 and y0i 1(0) = y0 0i 1(0) = = y(s i+1) 0i 1 (0) = 0 for i 2, we choose x1i, y0i satisfying x01i = Ax1i By0i; Gx1i = (GB)y0i 1; B0x1i(0) = 0: By the same argument as before, we obtain that x1i satis es (2.71c) for 2 i s+1, and y0i(0) = y0 0i(0) = = y(s i) 0i (0) = 0; s J + 1: Also, By0i is smooth. Next we denote the asymptotic solution x1s+1 = x10 + x11 + + sx1s + s+1x1s+1 and z1s+1 = x1 x1s+1: Then z0 1s+1 + Pz1s+1 = Az1s+1 + s+2( x01s+1 +Ax1s+1); z1s+1(0) = 0 Let u1s+1 = Rz1s+1 and w1s+1 = Pz1s+1. Hence, we have (cf. (2.38)) z1s+1 = Su1s+1 + w1s+1 and u01s+1 = (RA+R0)Su1s+1 + (RA +R0)w1s+1 +O( s+1); Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations49 w0 1s+1 + w1s+1 = (PA+ P 0)Su1s+1 + (PA+ P 0)w1s+1 +O( s+2); u1s+1(0) = 0; w1s+1(0) = 0: Using Lemma 2.2, we get w1s+1 = O( s+2) and u1s+1 = O( s+1), i.e. z1s+1 = O( s+1): Therefore, x1 = x10 + x11 + + sx1s +O( s+1): (2.76) Noting x10 = xe, we thus obtain x1 xe = O( ): (2.77) Then, by using (2.35),(2.76),(2.71),(2.72a) and (2.74a), it follows that By1 = By0 + 1 B(GB) 1(Gx1 + r) By1 = By0 + 1 (Px10 +B(GB) 1r + Px11 + + sPx1s +O( s+1)) = Bye + By01+ + s 1By0s 1 +O( s) (2.78) or By1 Bye = O( ): (2.79) Now we look at the second iteration s = 2 of (2.34), (2.35): x02 +B(GB) 1(Gx2 + r) = Ax2 By1 + q; with initial conditions (2.70). Let x2 = x20 + x21 + 2x22 + : Noting that (2.78) gives us a series expansion for By1 we obtain, B(GB) 1Gx20 = B(GB) 1r; (2.80a) B(GB) 1Gx21 = x020 +Ax20 Bye + q; (2.80b) B(GB) 1Gx2i = x02i 1 +Ax2i 1 By0i 1; 2 i s+ 1 (2.80c) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations50 Again, (2.80a) satis es initial conditions (2.70) and (2.80b) and (2.80c) satisfy the corresponding homogeneous ones. As the case of s = 1, we choose x20 = xe. We thus have B(GB) 1Gx21 = 0: Then x21 is constructed to satisfy x021 = Ax21 By11; (2.81a) Gx21 = 0; (2.81b) B0x21(0) = 0: (2.81c) Obviously x21 = 0 since (2.81) is uniquely solvable for x21 by Lemma 2.1. In general, similarly to the case of s = 1, we choose x2i satisfying x02i = Ax2i By1i; (2.82a) Gx2i = (GB)(y0i 1 y1i 1); (2.82b) B0x2i(0) = 0: (2.82c) for 2 i s+ 1. By applying Lemma 2.2 and the same argument as in the case of s = 1 we get x2 = xe + x21 + 2x22 + + sx2s +O( s+1) (2.83) or x2 xe = O( 2): (2.84) Then, using (2.35),(2.80),(2.81a),(2.82a), (2.83) and (2.78), we conclude By2 = By1 + 1 B(GB) 1(Gx2 + r) = Bye + 2By12 + + s 1By1s 1 +O( s) (2.85) or By2 Bye = O( 2) (2.86) Chapter 2. Sequential Regularization Methods for Di erential Algebraic Equations51 We can repeat this procedure, and, by induction, complete the proof for s J+1. 2 Chapter 3 SRM for Nonlinear Problems In the previous chapter we derived the sequential regularization method and gave a detailed continuous and discrete analysis for linear index-2 DAEs. The method relates to a combination of stabilization and penalty-like methods. In this chapter we extend the method to nonlinear index-2 and index-3 problems ( = 1 and = 2 in (2.15)), including constrained multibody systems. A number of variants are proposed, and particularly e ective methods are singled out in certain circumstances. All results obtained here are certainly applicable to the linear case. The chapter is organized as follows: In x3.1 we consider problems without constraint singularities. Two SRM variants are discussed. One variant involving dg dt (corresponding to (2.32), (2.33)) leads to nonsti problems. Taking E = I is particularly attractive. The other variant, corresponding to (2.32), (2.33) with 1 = 0, does not involve dg dt . The choice E = I, if possible (otherwise one can choose E = (GB)T ), makes the computation particularly simple. Problems with constraint singularities are considered in x3.2. The SRM corresponding to (2.34), (2.35) is proposed for such problems. This variant works well in practice, but our proofs to date extend only to the linear case. In x3.3 we analyze and discuss various methods for index-3 problems. A number of SRM variants are possible, combining regularization with Baumgarte's stabilization or with invariant stabilization [8]. Of particular interest, in case of no constraint singularity, are the methods (3.47) and (3.33){(3.35) which corresponds to invariant 52 Chapter 3. SRM for Nonlinear Problems 53 stabilization. The choice E = I leads to particularly simple iterations. A corresponding convergence result is given in Theorem 3.3. In case of a possible constraint singularity, the SRM (3.41) is recommended. These methods are reformulated in x3.4 for the special case of multibody systems with holonomic constraints. The \winning" methods are (3.44)-(3.45) with E = I for the nonsingular case and (3.46)-(3.47) for the case where the constraint Jacobian may have isolated rank de ciencies. In x3.5 we report the results of numerical experiments con rming our theoretical predictions and demonstrating the e ectiveness of the proposed methods. 3.1 Nonlinear, Nonsingular Index-2 Problems The nonlinear index-2 DAE ( = 1 in (2.15)) reads: x0 = f(x; t) B(x; t)y; (3.1a) 0 = g(x; t); (3.1b) where f , B and g are su ciently smooth functions of (x; t) 2 Rnx [0; tf ], and y 2 Rny . We consider this DAE subject to nx ny boundary conditions b(x(0); x(tf)) = : (3.2) These boundary conditions are assumed to yield a unique1 and bounded solution for the ODE (3.1a) on the manifold given by (3.1b). Concretely, if we were to replace (3.1b) by its di erentiated form (denoting G = @g @x) 0 = Gx0 + gt (= dg dt ) (3.3a) g(x(0); 0) = 0 (3.3b) 1locally unique, or isolated solution in a su ciently large neighborhood would su ce. Chapter 3. SRM for Nonlinear Problems 54 and use (3.3a) in (3.1a) to eliminate y and obtain nx ODEs for x, then the boundary value problem for x with (3.2) and (3.3b) speci ed has a unique solution. In the initial value case (i.e., when b is independent of x(tf)), this means that (3.2) and (3.3b) can be solved uniquely for x(0). In this section, we consider the case where GB is nonsingular. Generalizing the idea in x2.4, we have the following SRM formulation for the nonlinear index-2 DAEs (3.1), (3.2): for s = 1; 2; : : : ; x0s = f(xs; t) B(xs; t)ys; (3.4) where ys = ys 1 + 1 E(xs; t)( 1 d dtg(xs; t) + 2g(xs; t)); (3.5) subject to the boundary conditions (3.2) and (3.3b). Note that y0(t) is a given initial iterate which we assume is su ciently smooth and bounded and that > 0 is the regularization parameter. The regularization matrix E is nonsingular and has a uniformly bounded condition number; possible choices are E = I, E = (GB) 1 and others (e.g. E = (GB)T , cf. [11, 94]). We note that if we take y0 y then x1 x, where x and y are the solution of (3.1). If we take y0 0, then one SRM iteration is the usual penalty method (cf. [84, 91, 69]). As customary for the penalty method, we assume: Assumption 3.1 The problem (3.4), (3.5),(3.2),(3.3b) has a unique solution and the solution is bounded if ys 1 is bounded. Assumption 3.1 is generally true for initial value problems. For general boundary value problems, we expect that it would hold for most practical cases since (3.4) (with (3.5) plugged in) may be seen as a perturbed problem of (3.1) according to the proof of Theorem 3.1 (see below), where the perturbation and its rst derivative are both small if is small. Chapter 3. SRM for Nonlinear Problems 55 To analyze the SRM, we assume the following perturbation inequality: For 0 t tf , kx̂(t) x(t)k M max 0 tf(j ( )j+ j 0( )j); (3.6a) kŷ(t) y(t)k M max 0 tf(j ( )j+ j 0( )j); (3.6b) where k k is some lp norm (say, the maximum norm), and x̂ and ŷ satisfy the following perturbed version of (3.1): x̂0 = f(x̂; t) B(x̂; t)ŷ; (3.7a) 0 = g(x̂; t) + (t) (3.7b) with the same boundary conditions as (3.2). For initial value problems, (3.6) has been proved in [58], pp. 478-481. It is actually the de nition of the perturbation index introduced in [58]. Furthermore, (3.6) also holds for boundary value problems if we impose some boundedness conditions on the corresponding Green's function (cf. [14]). The case 1 6= 0 in (3.5) is su ciently di erent from the case 1 = 0 to warrant a separate treatment. 3.1.1 The case 1 = 1 Now we estimate the error of the sequential regularization method (3.4)-(3.5). We prove a theorem which says that the error after s SRM iterations is O( s) (i.e., each iteration improves the error by O( )) everywhere in t. This result coincides with that of the linear case. Theorem 3.1 Let all functions in the DAE (3.1) be su ciently smooth and the above assumptions hold. Then, for the solution of iteration (3.4),(3.5) with 1 6= 0, we have Chapter 3. SRM for Nonlinear Problems 56 the following error estimates: xs(t) xe(t) = O( s); ys(t) ye(t) = O( s); for 0 t tf and s 1. Proof: Let vs = g(xs; t). Then, from (3.4), v0 s = G(xs; t)x0s + gt(xs; t) = G(xs; t)f(xs; t) G(xs; t)B(xs; t)ys + gt(xs; t): Using (3.5), we thus have ( (GBE) 1 + I)v0 s + 2vs = (GBE) 1(Gf + gt) E 1ys 1; (3.8a) vs(0) = 0: (3.8b) Therefore it is not di cult to get vs = g(xs; t) = O( ); v0 s = g(xs; t)0 = O( ); (3.9) if ys 1 is bounded (which implies that xs is bounded from Assumption 3.1). For s = 1, we havex01 = f(x1; t) B(x1; t)y1 g(x1; t) = O( ); g(x1; t)0 = O( ) since y0 is chosen to be bounded. From assumption (3.6), we immediately get x1 xe = O( ); y1 ye = O( ): (3.10) Then it is easy to see that y1 is bounded. So for s = 2, we obtain x02 = f(x2; t) B(x2; t)y2 g(x2; t) = O( ); g0(x2; t) = O( ) Chapter 3. SRM for Nonlinear Problems 57 By using assumption (3.6) again, this yields x2 xe = O( ): Hence it can be veri ed, by substituting (3.3a),(3.1a) for the exact solution, that the right hand side of (3.8a) becomes O( 2). So, from (3.8), we can get g(x2; t) = O( 2); g0(x2; t) = O( 2): Applying assumption (3.6), it follows that x2 xe = O( 2); y2 ye = O( 2): (3.11) This also gives the boundedness of y2. We can repeat this procedure, and, by induction, conclude the results of the theorem. 2 From (3.8) it is clear that there is no sti ness here, so we can choose > 0 very small, so small in fact that one SRM iteration would su ce for any desired accuracy, and discretize the regularized ODE, possibly using a nonsti method like explicit Runge-Kutta. This gives a modi ed penalty method [I + 1BEG]x0 = f By0 1BE(gt + 2g) (3.12) where B;E; g etc, all depend on x, with the subscript s = 1 suppressed. For the choice E = (GB) 1, let P = BEG = B(GB) 1G be the associated projection matrix. Multiplying (3.12) by 1 1+ 1P and by I P , respectively, and then adding together, we have x0 = f 1 1 + 1By0 1 1 + 1B(GB) 1[Gf + gt + 2g] Thus the iteration obtained is similar to Baumgarte's stabilization x0 = f B(GB) 1[Gf + gt + 2g] (3.13) Chapter 3. SRM for Nonlinear Problems 58 In fact, the single SRM iteration tends to (3.13) in this case when ! 0. Indeed, the parameter 2 is the usual Baumgarte parameter, and choosing 2 > 0 obviously makes equation (3.8a) asymptotically stable for the drift vs. As indicated in [12], for both of these methods we can apply post-stabilization instead, i.e. take 2 = 0 but stabilize after each discretization step [8, 9]. For reasons of computational expense, it may be better to choose E = I in (3.12). The iteration obtained is simple, although a possibly large matrix (with a special structure) must be \inverted". Example 3.1 The choice of E = I is utilized in Chapter 4 (see also [80]) for the time-dependent, incompressible Navier-Stokes equations governing uid ow. The advantage gained is that no treatment of pressure boundary conditions is needed, unlike methods based on Baumgarte-type stabilizations which lead to the pressure-Poisson equation. 2 3.1.2 The case 1 = 0 For this case the drift equation (3.8) is clearly sti for 0 < 1. As in x2.5, we denote J such that y0(0) = ye(0); y0 0(0) = y0 e(0); : : : ; y(J) 0 (0) = y(J) e (0); (3.14) where J = 1 if y0(0) 6= ye(0), then we can prove the same result as Theorem 3.1 for s J + 1. Note that we may choose y0 satisfying (3.14) for some m 0 by expressing y in terms of x at t = 0 for initial value problems. But this starting procedure generally does not work for boundary value problems. Hence we state and prove the theorem for initial value problems and comment on the boundary value case following the proof. Chapter 3. SRM for Nonlinear Problems 59 Theorem 3.2 Let the assumptions of Theorem 3.1 plus (3.14) hold. In addition, suppose that the matrix function E(x; t) has been chosen so that GBE is positive de nite. Then, for the solution of iteration (3.4),(3.5) with 1 = 0, we have the following error estimates: xs(t) xe(t) = O( s); ys(t) ye(t) = O( s); for 1 s J + 1 and 0 t tf . Proof: We derive the result for the case s J + 1 = 2. Following the proof, we will comment on additional generalizations. The key is again the basic drift equation (3.8), which we rewrite here as v0 s + (GBE)vs = (Gf + gt GBys 1); (3.15a) vs(0) = 0: (3.15b) where quantities are evaluated as before, at (xs; t), unless otherwise noted. For s = 1, given the boundedness of y0 we obtain as before v1 = O( ) To obtain a similar result for v0 1, however, a di erent procedure from that of Theorem 3.1 is needed. Note that at t = 0, the condition y0(0) = y(0) implies (Gf + gt GBy0)jt=0 = 0 Hence from (3.15a), v0 1(0) = 0. Di erentiating (3.15a) with respect to t and using v1 = O( ), we get v00 1 + (GBE)v0 1 = O( ); v0 1(0) = 0 and this yields v0 1 = O( ) Chapter 3. SRM for Nonlinear Problems 60 From assumption (3.6) we then get (3.10). Subtracting (3.1a) from (3.4) and using (3.10) gives also x01 = x0 +O( ) and boundedness of y0 1 is obtained from a di erentiation of (3.5). For s = 2, given the boundedness of y1 (by (3.10)) and y1(0) = y(0), we get as for s = 1 v2 = O( ); v0 2 = O( ) and hence also x2 = x+O( ) This yields that the right hand side of (3.15a) is O( 2), so v2 = O( 2) Now comes the delicate part. To obtain an O( 2) estimate also for v0 2, so that the estimate (3.6) can be used to complete the proof, we di erentiate the drift equation again, obtaining v00 2 + (GBE)v0 2 = O( 2) + (Gf + gt GBy1)0 and v0 2(0) = 0 obtained as for the s = 1 case. We are then left to show that F (x2; t) := (Gf + gt GBy1)0 = O( ) (3.16) For this purpose we must estimate v00 1 rst. Using the condition y0 0(0) = y0(0), and also x1(0) = x(0), x01(0) = x0(0) (obtained from (3.4)), we can obtain (Gf(x1; t) + gt(x1; t) GB(x1; t)y0)0jt=0 = 0 Hence v00 1(0) = 0 from (3.15a) once di erentiated. Di erentiating (3.15a) twice we now obtain precisely as when estimating v0 1 above, v00 1 = O( ) Chapter 3. SRM for Nonlinear Problems 61 The boundedness of all needed quantities can also be obtained in the same way as before. Finally, we note v0 1 = Gx01 + gt(x1; t) = Gf(x1; t) + gt(x1; t) GB(x1; t)y1 Ready to show (3.16), we now write F (x2; t) = [(Gf(x2; t) Gf(x1; t))+(gt(x2; t) gt(x1; t))+(GB(x1; t) GB(x2; t))y1 v0 1]0 Our previous estimates allow the conclusion that x2 = x1 + O( ), x02 = x01 + O( ), hence we can nally conclude the estimate (3.16) and obtain the result of the theorem for s = 2. The proof proceeds in a similar manner for larger J . Generally, one needs the estimate v(j) 1 = O( ), 1 j J + 1, and this necessitates (3.14). 2 Remark 3.1 The convergence result holds for all s (i.e. also for s > J+1, assuming su cient smoothness) away from an initial layer of size O( ) in t. This is so because E is chosen so that we can express the solution for small as a smooth outer solution which is bounded in terms of the right hand side as before, plus an initial layer of width O( ). Conditions (3.14) then ensure that the layer error is bounded by O( J+1) for the rst J + 1 iterations. 2 Remark 3.2 For boundary value problems, there is no obvious technique to ensure J > 1. For a given J , the results of Theorem 3.2 and Remark 3.1 can be extended as in Chapter 2. This requires a di erent proof technique, though. Basically, an asymptotic expansion for xs and ys is constructed, where the rst term is the exact solution x; y. This latter proof technique follows more along traditional singular perturbation lines (see [112, 89, 69]), and is not as close to Theorem 3.1 and to DAE concepts. 2 Chapter 3. SRM for Nonlinear Problems 62 Taking 2 = 1 without loss of generality, we obtain the iteration x0s = f Bys 1 1BEg(xs; t) (3.17) This is a singular, singularly perturbed problem (so should not be taken extremely small compared to machine precision even if a sti solver is being used). If GB is positive de nite then we may choose E = I, and this yields a very simple iteration in (3.17) which avoids the inversion necessary in stabilization methods like Baumgarte's. However, if an explicit discretization method of order p is contemplated then approximately p SRM iterations like (3.17) are needed, because one must choose = O(h), where h is the step size. 3.2 Nonlinear, Singular Index-2 Problems In this section we consider the nonlinear index-2 problem (3.1) with an isolated singular point t?, i.e. GB is singular at t?. For simplicity, we assume that B and g are independent of t. Denote P (x) = B(GB) 1G. Motivated by constrained multibody systems (see Example 2.1), we assume P (x) to be di erentiable in t, but @P @x (x) may be unbounded. For this reason, we consider only the case 1 = 0 in this section (cf. Chapter 2 for the linear case). In the drift equation (3.8) we then have essentially the singularly perturbed operator v0 + GBEv to consider. The choices of E = I or E = (GB)T yield a turning point problem (i.e., at least one of eigenvalues of the matrix GBE vanishes at the point t?), which complicates the analysis, even in the linear case , and degrades the numerical performance as well in our experience. Therefore, we choose E = (GB) 1. In the sequel we will be careful to evaluate the e ect of E only when its singularity limit is well-de ned, as e.g. in P (x). A direct generalization of the linear case in Chapter 2 would give the SRM formulation (3.4), (3.5) where instead of updating y (because y may be unbounded at Chapter 3. SRM for Nonlinear Problems 63 t ) we update By by B(xs)ys = B(xs 1)ys 1 + 1 B(xs)(G(xs)B(xs)) 1g(xs): (3.18) However, (3.18) needs to be modi ed, since we may haveRangeB(xs) 6=RangeB(xs 1). So we use the projection P (xs) to move from RangeB(xs 1) to RangeB(xs). Then we consider the following SRM formulation for singular problems: x0s = f(xs; t) B(xs)ys; (3.19a) B(xs)ys = P (xs)B(xs 1)ys 1 + 1 B(xs)(G(xs)B(xs)) 1g(xs); (3.19b) where xs satis es the boundary condition (3.2). If the assumptions given at the beginning of x3.1 and in Theorem 3.2 remain valid, then the result of Theorem 3.2 still holds. Unfortunately, for the singular problem, assumption (3.6) may not be true in general. To see this, consider one iteration, i.e. s = 1. The accuracy for the approximation of x depends on the extent that the bound (3.6a) holds. Numerical experiments show that we can get a good approximation of x near the singularity. But the situation for By is worse, and the bound (3.6b) often does not hold. Indeed, assume for the moment that we have a good, smooth approximation of x, say xs = x̂, i.e. (3.7) holds with ; 0 = O( ), and B(x̂)ŷ is de ned by (3.19b) for some B(xs 1)ys 1. From (3.7) we have B(x̂)ŷ = P (x̂)f(x̂; t) + ; (3.20) where = B(x̂)(G(x̂)B(x̂)) 1 0. It is not di cult to nd that the exact B(x)y from (3.1) satis es B(x)y = P (x)f(x; t): (3.21) Yet, even if is small, B(x̂)ŷ may not be a good approximation of By because @P @x may be unbounded at the singular point so that P (x̂) is not a good approximation of P (x). Chapter 3. SRM for Nonlinear Problems 64 Example 3.2 In (3.1) let x = (x1; x2), g(x) = cosx1 cos x2, and G = BT = ( sinx1 sinx2 ). Then P (x) = (sin2 x1 + sin2 x2) 1 0@ sin2 x1 sinx1 sinx2 sinx1 sin x2 sin2 x2 1A. Clearly, at a singular point x = (0; ), the value of P depends on the direction from which it is approached. Thus, @P @x is unbounded, even though P is a di erentiable function of t. Further letting f = (sinx2 sin x1) 1 ( sinx2 2 sin x2 sin x1 )T , and given the initial conditions x1(0) = =2; x2(0) = =2, the exact solution is x(t) = ( t =2 t+ =2 )T ; y = (sin x2 sinx1) 1 Thus, as t crosses t = =2, y(t) becomes unbounded, but By = (sinx2 sinx1) 1 ( sinx1 sinx2 )T remains bounded. However, it is easy to perturb x(t) slightly and smoothly in such a way that the perturbed By becomes unbounded near t = t , still satisfying (3.7) with a small . 2 Note that for the linear model problem (see Chapter 2), P P (t) is independent of x. Hence we do not have the above di culty in the linear case. For the nonlinear problem, the accuracy near the singular point is reduced and it no longer behaves like O( s) for more than one iteration. However, we do expect O( s) accuracy away from the singular point, assuming that no bifurcation or impasse point is encountered by the approximate solution . 3.3 SRM for Nonlinear Higher-index Problems We now generalize the SRM to the more general problem (2.15). In particular, we consider the index-3 problem ( = 2). The Euler-Lagrange equations for multibody systems with holonomic constraints yield a practical instance of the problem. The Chapter 3. SRM for Nonlinear Problems 65 SRM formulations presented in this section are easy to generalize for more general problems (2.15) (index + 1). The index-3 problem reads: x00 = f(x; x0; t) B(x; t)y; (3.22a) 0 = g(x; t); (3.22b) with given 2(nx ny) boundary conditions, b(x(0); x(tf); x0(0); x0(tf)) = 0: (3.23) The meaning of G, B and the stabilization matrix E below remain the same as in the index-2 problems considered in previous sections. 3.3.1 The case of nonsingular GB We rst use the idea from previous sections, viz. a combination of Baumgarte's stabilization with a modi ed penalty method, to derive the SRM for the nonlinear index-3 problem (3.22). Then we apply a better stabilization [8] to generate a new SRM which is expected to have better constraint stability. Finally, we seek variants which avoid evaluation of complicated terms in the second derivative of the constraints. First consider, instead of (3.22b), the Baumgarte's stabilization 1 d2 dt2g(x; t) + 2 d dtg(x; t) + 3g(x; t) = 0; (3.24a) g(x(0); 0) = 0; d dtg(x(0); 0) = 0; (3.24b) where j; j = 1; 2; 3 are chosen so that the roots of the polynomial ( ) = 3 X j=1 j 3 j are all negative. Following the procedure of previous sections, we can write down an SRM for (3.22): for s = 1; 2; : : : and y0 given, x00 s = f(xs; x0s; t) B(xs; t)ys; (3.25) Chapter 3. SRM for Nonlinear Problems 66 where xs satis es boundary conditions (3.23) and (3.24b) and ys is given by ys = ys 1 + 1 E(xs; t)( 1 d2 dt2g(xs; t) + 2 d dtg(xs; t) + 3g(xs; t)): (3.26) It is not di cult to repeat the approach of x3.1 for the present case. Under assumptions similar to the index-2 case, i.e. (3.6) with a change to include 00( ) at the right hand side (cf. [58]) and Assumption 3.1 with the addition that the derivative of the solution is also bounded, we readily obtain extensions of Theorems 3.1 and 3.2 for the cases 1 6= 0 and 1 = 0 (with 2 6= 0), respectively. We do not allow 1 = 2 = 0 since in this case equations (3.25),(3.26) have di erent asymptotic properties. Note that the SRM (3.25),(3.26) with 1 = 0 avoids computing gxx; however, the iteration obtained now calls for solving problems which become sti when gets small, and to avoid gxx one should use a non-sti discretization method. This formulation with E = I and 1 6= 0 is the same as that proposed in [20, 19] using the augmented Lagrangian method. Another way to generalize the SRM to higher index problems is based on invariant stabilization. Its advantages over Baumgarte's stabilization have been discussed in [8, 9]. We thus prefer the way based on this stabilization. Theoretical evidence is also mentioned in Remark 3.4. We rst describe this new stabilization. By two direct di erentiations of the constraints (3.22b), we can eliminate y and get an ODE x00 = ~ f(x; x0; t); (3.27) for which the original constraint (3.22b) together with its rst derivative give an invariant. The idea of the method is to reformulate the higher index DAE (3.22) as a rst order ODE (cf. (3.27)): z0 = f̂ (z; t) (3.28) with an invariant 0 = h(z; t); (3.29) Chapter 3. SRM for Nonlinear Problems 67 where z = 0@ z1 z2 1A = 0@ x x0 1A ; f̂(z; t) = 0@ z2 ~ f (z; t) 1A ; h(z(t); t) = 0@ g(x(t); t) d dtg(x(t); t) 1A (3.30) and to consider the stabilization families z0 = f̂(z; t) F (z; t)h(z; t); (3.31) where F = D ~ E for some appropriate matrix functions D and ~ E such that ~ E and HD are nonsingular and H = hz. The ODE (3.31) coincides with Baumgarte's stabilization for the index-2 problem (3.1) with D = B and ~ E = E = (HD) 1. One choice for D is D = HT , but others will be mentioned below. Note that (3.31) has the same solution as the original problem (3.22) for any parameter value . Although the method has better constraint stabilization, both the evaluation of ~ f and that of H involve gxx which may be complicated to calculate in practice. Next, we present an SRM method based on invariant stabilization which avoids the computation of ~ f . In fact, we can avoid gxx altogether using the new stabilization. If we do not eliminate y by di erentiations, f̂ (z; t) in the stabilization (3.31) becomes f̂(z; t) = 0@ z2 f(z; t) B(z1; t)y 1A : (3.32) Since y is not known in advance, we use an iterative SRM procedure to calculate y as in [20, 11]. The solutions of the iterative procedure no longer satisfy (3.22) precisely. Hence the iterative procedure has to be a regularization procedure and the parameter in (3.31) is changed to = 1 to emphasize that it must be chosen su ciently large. These lead to the following SRM formulation (for simplicity of notation, we only consider the special case where B and g are independent of t): z0 s = 0@ z1s z2s 1A0 = 0@ z2s f(zs; t) B(z1s)ys 1 1A 1 F (zs)h(zs); (3.33) Chapter 3. SRM for Nonlinear Problems 68 where zs satis es boundary conditions (3.23), (3.24b) and h = (g(z1); G(z1)z2)T . Thus the Jacobian of h is H = 0@ G(z1) 0 L(z) G(z1) 1A ; where L = zT 2 gxx(z1): We choose D and ~ E so that F = BE 0@ I 0 0 I 1A = 0@ BE 0 0 BE 1A (3.34) where, as in x3.1, E is chosen such that GBE has non-negative eigenvalues. Updating y by ys = ys 1 + 1 E(z1s)G(z1s)z2s (3.35) implies that the second part of the original index-3 system holds exactly, i.e. z0 2s = f(zs; t) B(z1s)ys: Next we analyze the convergence of (3.33){(3.35). Again we assume that the solutions of (3.33), (3.23), (3.24b) exist uniquely and are bounded if ys 1 is bounded (see Assumption 3.1 ). Assumption (3.6) changes slightly: We rst rewrite the system (3.22) as z0 1 = z2; (3.36a) z0 2 = f(z; t) B(z1)y; (3.36b) 0 = g(z1): (3.36c) Then we assume the following perturbation bound, kẑ(t) z(t)k M max 0 tf(j ( )j+ j 0( )j+ j 00( )j+ j ( )j+ j 0( )j); (3.37a) kŷ(t) y(t)k M max 0 tf(j ( )j+ j 0( )j+ j 00( )j+ j ( )j+ j 0( )j); (3.37b) Chapter 3. SRM for Nonlinear Problems 69 where ẑ and ŷ satisfy a perturbed problem of (3.36), ẑ0 1 = ẑ2 + (t); (3.38a) ẑ0 2 = f(ẑ; t) B(ẑ1)ŷ; (3.38b) 0 = g(ẑ1) + (t); (3.38c) with the same boundary conditions (3.23). Again, for initial value problems, (3.37) can be easily proved by following the technique presented in [58], and this can be extended for boundary value problems as well. Similarly to the proof of Theorem 3.1, let h(zs) = 0@ vs ws 1A, where vs = g(z1s), ws = G(z1s)z2s. From (3.33), we get the drift equations (cf. (3.15)) v0 s = GBE(zs)vs + ws (3.39a) v00 s = GBE(zs)v0 s LBE(zs)vs + [Lz2s +Gf(zs) GB(zs)ys 1] (3.39b) with the initial conditions vs(0) = 0, ws(0) = 0. Applying (3.39a) and then (3.39b) for s = 1, we obtain v1 = O( ); w1 = O( ). Then (3.39a) further yields v1 = O( 2); v0 1 = O( ) Comparing (3.38) with (3.33){(3.35), we have to bound = v1; = 1BEv1 and their derivatives appearing in (3.37). We already have that = O( 2); 0 = O( ); = O( ) The procedure that follows continues to be similar to the one employed in the proof of Theorem 3.2, so we only sketch it here for s = 1. From (3.39a) we obtain v0 1(0) = 0. Using the condition y0(0) = y(0) gives [Lz21+Gf(z1) GB(z1)y0]jt=0 = 0 Chapter 3. SRM for Nonlinear Problems 70 so, from (3.39b), also w0 1(0) = 0. Di erentiating (3.39b) we get w00 1 +GBE(z1)w0 1 = O( ); w0 1(0) = 0 hence w0 1 = O( ). Di erentiating (3.39a) we next have v00 1 +GBE(z1)v0 1 = O( 2); v0 1(0) = 0 so v0 1 = O( 2). This implies 0 = O( ); 00 = O( ) We can now use (3.37) and obtain the desired conclusion for s = 1, z1 = z +O( ); y1 = y +O( ); where fz; yg is the exact solution of the index-3 problem. Then, continuing to follow the proof procedure of Theorem 3.2, we obtain: Theorem 3.3 Let all functions in the DAE (3.22) be su ciently smooth and assume the above assumptions (particularly (3.37)) hold. Assume in addition that y0 satis es (3.14). Then, for the solution of iteration (3.33){(3.35), the following error estimates hold: zs(t) ze(t) = O( s); (3.40a) ys(t) ye(t) = O( s) (3.40b) for 1 s J + 1. Remark 3.3 Extensions of this theorem to the boundary value case and to s > J+1 away from an initial layer are possible, similarly to the extensions for Theorem 3.2 contained in Remarks 3.1 and 3.2. 2 Chapter 3. SRM for Nonlinear Problems 71 Remark 3.4 We note that, unlike Proposition 2.2 of [8], we do not assume kH(z)f̂(z)k2 0kh(z)k2 to discuss the stability and accuracy of the constraints. Also, from (3.39), we see the di erence of the constraint stability or accuracy between SRM formulations based on Baumgarte's stabilization and the new stabilization. For the former, we only have v0 1 = G(z11)z0 11 = G(z11)z21 = w1 So if we obtain w1 = O( ) then v1 = O(t ). This can be much worse than what we get from (3.39a). 2 3.3.2 The case for constraint singularities For the singular case GB may be singular at some isolated point t as described in the previous sections. The situation here is similar to that for index-2 problems. An examination of the drift equations (3.39) suggests that here, too, the choice E = (GB) 1 is preferable to E = I or E = (GB)T . The iteration for ys is modi ed as well. Still assuming for simplicity that g and B do not depend explicitly on t, this gives, in place of (3.33){(3.35) the iteration z0 1s = z2s 1 B(GB) 1g(zs) (3.41a) z0 2s = f(zs; t) ŷs (3.41b) ŷs = P (zs)ŷs 1 + 1 P (zs)z2s (3.41c) Also, as indicated in x3.2 for index-2 problems, we cannot expect O( s) approximation near the singular point any more. But we do expect that (3.40) holds away from the singular point, because the singularity is in the constraint and the drift manifold is asymptotically stable (following our stabilization). A numerical example Chapter 3. SRM for Nonlinear Problems 72 in x3.5 will show that we do get improved results by using SRM iterations for the singular problem. 3.4 SRM for Constrained Multibody Systems Constrained multibody systems provide an important family of applications of the form (3.22) and (3.1). We consider the system q0 = v (3.42a) M(q)v0 = f(q; v) G(q)T (3.42b) 0 = g(q) (3.42c) where q and v are the vectors of generalized coordinates and velocities, respectively; M is the mass matrix which is symmetric positive de nite; f(q; v) is the vector of external forces (other than constraint forces); g(q) is the vector of (holonomic) constraints; is the vector of Lagrange multipliers; and G(q) = d dqg. For notational simplicity, we have suppressed any explicit dependence of M , f or g on the time t. We rst consider the problem without singularities. Corresponding to (3.22) in x3.3, we have B = M 1GT , so GB = GM 1GT . Other quantities like h and H retain their meaning from the previous section. In some applications it is particularly important to avoid terms involving gxx, since its computation is somewhat complicated and may also easily result in mistakes and rugged terms. So [9] suggests post-stabilization using the stabilization matrix F =M 1GT (GM 1GT ) 10@ I 0 0 I 1A (3.43) twice, instead of involving H, at the end of each time step or as needed. They nd that this F performs very well in many applications. However, while this stabilization avoids the gxx term in F , gxx is still involved in obtaining ~ f , although only through Chapter 3. SRM for Nonlinear Problems 73 matrix-vector multiplications (see (3.27)). The SRM formulation (3.33){(3.35) enables us to avoid the computation of ~ f in the absence of constraint singularities. For the multibody system (3.42) we write the iteration as follows: For s = 1; 2; : : :, nd fqs; vsg by q0 s = vs 1 BE(qs)g(qs) (3.44a) v0 s = M 1f(qs; vs) B(qs) s 1 1 BEG(qs)vs (3.44b) Then update by s = s 1 + 1 EG(qs)vs: (3.45) It is easy to see that in this SRM formulation the gxx term is avoided completely. Moreover, since GM 1GT is positive de nite, we can choose E = I in (3.44),(3.45), obtaining a method for which Theorem 3.3 applies, which avoids computing (GM 1GT ) 1. Although it requires an iterative procedure, a small number of iterations (p if an explicit discretization method of order p is used) typically provide su cient accuracy. Numerical experiments will show the O( s) error estimate. Next we consider the singular problem, i.e. with the matrix GM 1GT being singular at some isolated point t , 0 < t < tf . A typical example of singular multibody systems is the two-link slider-crank problem (see Example 2.1 and Figure 2.1) consisting of two linked bars of equal length, with one end of one bar xed at the origin, allowing only rotational motion in the plane, with the other end of the other bar sliding along the x-axis. Various formulations of the equations of motion for this problem appear, e.g., in [60, 19, 11, 12, 94]. In our calculations we have used the formulation in Example 2.1 (see [11]), to make sure that the problem is not accidentally too easy. It consists of 6 ODEs and 5 constraints, with the last row of the Jacobian matrix G vanishing when the mechanism moves left through the point where both bars are upright ( 1 = 2 ; 2 = 3 2 , where xi; yi; i are the coordinates of the center of mass of the ith bar). Chapter 3. SRM for Nonlinear Problems 74 The last row of G vanishes at this point and a singularity is obtained. We note that the solution is smooth in the passage through the singularity with a nonzero velocity. When we attempt to integrate this system using a stabilization method like [8] which ignores the singularity, the results are unpredictable, depending on how close to the singular time point the integration process gets when attempting to cross it. In fact, radically di erent results may be obtained upon changing the value of the error tolerance. (Similar observations are made in [94].) In some instances a general purpose ODE code would simply be unable to \penetrate the singularity" and yield a solution which, after hovering around the upright (singular) position for a while, turns back towards the initial position (solid line in Figure 2.1). Such a pattern of motion may well look deceptively plausible. Methods which do not impose the constraints on the position level (e.g. methods consisting of di erentiating the constraints once and solving the resulting index-2 problem numerically, or of projecting only on the velocity-level constraint manifold) perform particularly poorly (cf. numerical results in [94]). This is easy to explain: The position-level constraint corresponds to ensuring that the two bars have equal length. If this is not strictly imposed in the process of numerical solution, inevitable numerical errors due to discretization may yield a model where the lengths are not close enough to being equal, and this leads to the lock-up phenomena described e.g. in [60], which have a vastly di erent solution pro le. We now wish to generalize the SRM to problem (3.42) with singularities since we have seen its success for the linear index-2 case. From the two-link slider crank problem, we nd that, although GM 1GT is singular at t , P (q) M 1GT (GM 1GT ) 1G and M 1GT (GM 1GT ) 1g are smooth functions of t for the exact solution or functions q satisfying the constraints, whileM 1GT (GM 1GT ) 1,M 1GT (GM 1GT ) 1Gq and the derivative dP (q) dq are not. Also, as indicated in [11], is no longer smooth, while B is since we assume the solution q to be su ciently smooth. We only include Chapter 3. SRM for Nonlinear Problems 75 terms which are (most possibly) smooth in the SRM formulation. Applying (3.41), we obtain the method q0 s = vs 1 M 1GT (GM 1GT ) 1g(qs); (3.46a) v0 s = M 1f(qs; vs; t) ̂s; (3.46b) ̂s = P (qs)̂s 1 + 1 P (qs)vs (3.47) As we indicated in x3.2, we do not expect O( s) accuracy near the singular point. However, we do expect that the SRM iteration would improve the accuracy and that we still expect to get O( s) accuracy away from the singular point. Numerical experiments in x3.5 will show such improvements. 3.5 Numerical Experiments We now present a few examples to demonstrate our claims in the previous sections. Throughout this section we use a constant step size h. To make life di cult we choose h when we can so that there is an i such that ti = t , namely, there is a mesh point hitting the singularity point t , for singular test problems. At a given time t, we use 0ex0 to denote the maximum over all components of the error in xs. Similarly, 0drift0 denotes the maximum residual in the algebraic equations. Example 3.3 Consider the DAE (2.16),(2.17) with f = 0@ 1 e t cos t+ et sin t1A ; B = 0@ x1 x21A g = 1 2(x21 + x22 e 2t sin2 t): subject to x1(0) = 1; x2(0) = 0. The exact solution is xe = ( e t; sin t ), ye = et. This is a problem without singularities. Chapter 3. SRM for Nonlinear Problems 76 Using an explicit second order Runge-Kutta method with h = 0:001 we test various choices of E and 1 (always taking 2 = 1) of the SRM formulation in x3.1. We list the computational results in Table 3.1. Observe that, for 1 6= 0, the SRM works well methods iteration error at ! t=.1 t=.5 t=1.0 1 = 1 1e-8 1 ex .11e-7 .94e-7 .19e-6 E = I drift .79e-8 .56e-7 .14e-6 1 = 1 1e-8 1 ex .11e-7 .92e-7 .18e-6 E = (GB)T drift .78e-8 .53e-7 .14e-6 1 = 1 1e-8 1 ex .11e-7 .95e-7 .19e-6 E = (GB) 1 drift .80e-8 .58e-7 .15e-6 Baumgarte ex .45e-6 .16e-6 .35e-6 drift .40e-6 .70e-7 .29e-6 1 = 0 5e-3 1 ex .60e-2 .11e-1 .11e-1 E = I drift .54e-2 .80e-2 .13e-1 2 ex .11e-3 .26e-3 .22e-3 drift .96e-4 .20e-3 .27e-3 3 ex .32e-5 .65e-5 .46e-5 drift .29e-5 .47e-5 .54e-5 4 ex .26e-6 .23e-6 .28e-6 drift .13e-6 .51e-7 .12e-6 1 = 0 5e-3 1 ex .70e-2 .12e-1 .13e-1 E = (GB)T drift .64e-2 .13e-1 .15e-1 2 ex .22e-3 .65e-3 .31e-3 drift .20e-3 .49e-3 .29e-3 3 ex .11e-4 .16e-4 .69e-5 drift .10e-4 .10e-4 .52e-5 4 ex .85e-6 .91e-7 .29e-6 drift .75e-6 .77e-6 .14e-6 1 = 0 5e-3 1 ex .51e-2 .66e-2 .10e-1 E = (GB) 1 drift .46e-2 .49e-2 .12e-1 2 ex .35e-4 .11e-3 .21e-3 drift .30e-4 .79e-4 .24e-3 3 ex .86e-6 .23e-5 .47e-5 drift .77e-6 .17e-5 .53e-5 4 ex .26e-6 .18e-6 .26e-6 drift .26e-7 .31e-7 .13e-6 Table 3.1: Errors for Example 3.3 using the explicit second order Runge-Kutta scheme for various choices of E. Its error is as good as Baumgarte's method whose parameter corresponds to the 2 of the SRM. For 1 = 0, we see that the error improves at a rate of about O( ) for various choices of E, including E = I. (Observe the errors at t = 1; the error situation near t = :1 is di erent because of an initial layer.) Such an error Chapter 3. SRM for Nonlinear Problems 77 improvement continues until the accuracy of the second order explicit Runge-Kutta method, i.e. O(h2), is reached. 2 The next two examples are problems with singularities. In the index-2 case of the Baumgarte stabilization the worst term is B(GB) 1gt for the type of the singularities in this paper. To show what happens when the Baumgarte method does not work well, we choose nonautonomous problems (i.e. gt 6= 0) as index-2 singular examples. Example 3.4 Consider the nonlinear DAE (2.16) with f = 0@ 1 + (t 12)et 2t+ (t2 1 4)et1A ; B = 0@ x1 x21A g = 12(x21 + x22 (t 1 2)2 (t2 1 4)2) subject to the initial condition x1(0) = 12 ; x2(0) = 14. The exact solution is xe = ( t 12 ; t2 1 4 ), ye = et. A singularity is located at t = 12 . Using this example we test the SRM formulations of x3.2. We list the computational results in Table 3.2, where we take h = = 0:001 for the case of 1 = 0, and h = 0:001; = 10 10 for the case of 1 6= 0, and use the explicit second order Runge-Kutta scheme to easily see the iteration improvement (Ij stands for results of the jth iteration). From Table 3.2, we see the error's deterioration for the Baumgarte method and the SRM with 1 6= 0. The SRM with 1 = 0 performs better in the singular case. 2 Next we try an example in which y is unbounded at the singularity. Example 3.5 Consider the nonlinear DAE (3.1) with f = 0@ x1 + x2 sin(t) (1 + 2t) 0 1A ; B = 0@ 0 x11A g = x21 + x1(x2 sin(t) 1 + 2t); Chapter 3. SRM for Nonlinear Problems 78 methods error at ! t=.1 t = .3 t = .5 t= .7 t= 1.0 1 = 1 ex .39e-6 .13e-5 .12e-3 .14e-3 .76e-4 drift .24e-6 .16e-6 .10e-7 .39e-6 .75e-6 1 = 0 (I1) ex .46e-3 .32e-3 .43e-4 .49e-3 .20e-2 drift .24e-3 .89e-4 .18e-8 .20e-3 .22e-2 1 = 0 (I2) ex .81e-6 .11e-5 .41e-5 .29e-5 .68e-5 drift .24e-6 .30e-6 .15e-10 .13e-5 .76e-5 1 = 0 (I3) ex .23e-6 .26e-6 .34e-6 .29e-6 .29e-6 drift .90e-9 .11e-8 .78e-13 .35e-8 .18e-7 1 = 0 (I4) ex .23e-6 .26e-6 .36e-6 .27e-6 .29e-6 drift .47e-11 .33e-11 .10e-12 .29e-11 .28e-10 Baumgarte ex .43e-6 .45e-6 .34e-3 .39e-3 .21e-3 drift .24e-6 .16e-6 .61e-7 .24e-6 .75e-6 Table 3.2: Example 3.4 { bounded y and singularity at t = :5 subject to the initial condition x1(0) = 1; x2(0) = 0. The exact solution is xe = ( 1 2t; sin t ), ye = cos t=(1 2t). Taking the same parameters and using the same method as before, we get the results listed in Table 3.3. Clearly, the SRM with 1 = 0 performs well for this situation, while methods error at ! t=.1 t = .3 t = .5 t= .7 t= 1.0 SRM ( 1 = 0) ex .40e-6 .25e-6 .14e-6 .46e-7 .60e-7 (I3) drift .25e-8 .76e-9 .16e-15 .28e-9 .40e-9 Baumgarte ex .49e-7 .15e-6 .93e+1 NaN NaN drift .39e-7 .59e-7 .52e+13 NaN NaN Table 3.3: Example 3.5 { unbounded y and singularity at t = :5 Baumgarte method blows up upon hitting the singularity. 2 Our next example tests the formulation (3.33){(3.35) or (3.44){(3.45) for index-3 problems. Example 3.6 This example is made up from Example 2 in [9] (see Figure 3.1), which describes a two-link planar robotic system. We use the notation of (1.11). Let Chapter 3. SRM for Nonlinear Problems 79 θ 1 2 x y 1 1 2 2 θ2 1 l l (x , y ) (x , y ) Figure 3.1: Two-link planar robotic system q = ( 1; 2)T and M = 0@m1l2 1=3 +m2(l2 1 + l2 2=3 + l1l2c2) m2(l2 2=3 + l1l2c2=2) m2(l2 2=3 + l1l2c2=2) m2l2 2=3 1A ; where l1 = l2 = 1, m1 = m2 = 3 and c2 = cos 2. The constraint equation is g(q) = l1 sin 1 + l2 sin( 1 + 2) = 0: We choose the force term f = 0@ (l1 cos 1 + l2 cos( 1 + 2)) cos t 3 sin t l2 cos( 1 + 2) cos t+ (1 3 2c2) sin t 1A which yields the exact solution 1 = sin t, 2 = 2 sin t and = cos t. Because M is (symmetric) positive de nite and B =M 1GT we can take E = I in the SRM formula (3.44){(3.45). Again we use the secondorder explicit Runge-Kutta scheme, and set h = 0:001; = 0:005. The results are listed in Table 3.4, where eq and ev stand for maximum errors in q and v = q0, resp., and pdrift and vdrift stand for drifts at the Chapter 3. SRM for Nonlinear Problems 80 position and velocity level, resp. We see that the accuracy is improved signi cantly by the rst two iterations. The third iteration is unnecessary here, because the error is already dominated by the Runge-Kutta discretization error. Qualitatively similar results are obtained for E = (GB)T and E = (GB) 1. More interestingly, though, for E = I we neither form nor invert GM 1GT , so a particularly inexpensive iteration is obtained. methods iteration error at ! t=.1 t=.5 t=1.0 E = I 5e-3 1 eq .41e-4 .66e-3 .26e-2 ev .75e-2 .74e-2 .69e-2 pdrift .22e-4 .28e-4 .22e-4 vdrift .49e-2 .41e-2 .27e-2 2 eq .13e-6 .66e-6 .36e-6 ev .19e-5 .81e-6 .20e-4 pdrift .42e-9 .13e-7 .17e-6 vdrift .91e-7 .21e-5 .21e-4 3 eq .10e-6 .58e-6 .12e-5 ev .86e-6 .10e-5 .16e-5 pdrift .96e-11 .60e-9 .48e-8 vdrift .10e-8 .99e-7 .59e-6 Table 3.4: Errors for Example 3.6 using SRM (3.45)-(3.46) 2Next we solve for the dynamics of the slider-crank mechanism described in Example 2.1. this is a nonlinear index-3 DAE with isolated, \smooth" singularities. Example 3.7 We take h = = 0:0001 and use the explicit second order RungeKutta method again. Singularities are located at ( 1; 2) = ( 2 ; 3 2 ) (i.e., each time the periodic solution passes this point). Corresponding to the case shown in [94], we choose 1(0) = 7 4 and 01(0) = 0 and compute 1 = 1 3 2 ; 2 = 2 + 2 ; 0 1 and 0 2. Using the formulation (3.47), (3.46), we calculate until t = 70 without any di culty (see Figure 3.2). Chapter 3. SRM for Nonlinear Problems 81 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 0 10 20 30 40 50 60 70 "theta1" "theta2" "dtheta1dt" "dtheta2dt" Figure 3.2: Solution for slider-crank problem with singularities We also list the drift improvement as a function of the SRM iteration in Table 3.5. iteration number position drift at t=30 velocity drift at t = 30 1 .669e-8 .671e-4 2 .730e-11 .731e-7 Table 3.5: Drifts of the SRM for the slider-crank problem If we use the SRM formulations considered in xx3.3 and 3.4 for problems withour singularities, or one of the usual stabilization methods with strict tolerances, the results become wildly di erent from the correct solution after several periods. Next we calculate the acceleration of the slider end in the horizontal direction under the initial data 1(0) = 4 and 01(0) = 2p2. The same problem was discussed in [19]. The result shown in [19] is not perfect since the maximum and minimum values in each period appear to di er. Our result looks better (see Figure 3.3). 2 Chapter 3. SRM for Nonlinear Problems 82 0 1 2 3 4 5 6 7 8 9 10 −80 −60 −40 −20 0 20 40 60 80 Time(s) X a cc el er at io n of c ra nk e nd ( m /s ^2 ) Figure 3.3: Acceleration of slider end Chapter 4 SRM for the Nonstationary Incompressible Navier-Stokes Equations 4.1 DAE Methods for Navier-Stokes Equations While a signi cant body of knowledge about the theory and numerical methods for DAEs has been accumulated, not much has been extended to partial di erentialalgebraic equations (PDAEs). The incompressible Navier-Stokes equations form, in fact, an example of a PDAE: to recall, these equations read ut + (u grad)u = u gradp+ f ; (4.1a) divu = 0; (4.1b) uj@ = b ; ujt=0 = a; (4.1c) in a bounded twoor three-dimensional domain and the time interval 0 t T . Here u(x; t) represents the velocity of a viscous incompressible uid, p(x; t) the pressure, f the prescribed external force, a(x) the prescribed initial velocity, and b(t) the prescribed velocity boundary values. The system (4.1) can be seen as a partial di erential equation with constraint (4.1b) with respect to the time variable t. Comparing with the DAE form (3.1) p corresponds to y, the grad operator corresponds to the matrix B and the div operator corresponds to the Jacobian matrix G. It is easily veri ed that (4.1) has index-2 since the operator divgrad = (corresponding to the matrix GB) is invertible (under appropriate boundary conditions). Indeed, the pressure-Poisson reformulation of (4.1) (see, e.g., [55]) corresponds to a direct index reduction of the PDAE, i.e. a di erentiation of the constraint with respect to t followed by substitution of ut from the momentum equations. 83 Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 84 In this chapter we propose and analyze a sequential regularization method (SRM) for solving the incompressible Navier-Stokes equations. The method is de ned as follows: with p0(x; t) an initial guess, for s = 1; 2; ; solve the problem (us)t grad( 1(divus)t + 2divus) + (us grad)us = us gradps 1 + f ; (4.2a) usj@ = b;usjt=0 = a; (4.2b) ps = ps 1 1 ( 1(divus)t + 2divus): (4.2c) This method is an extension of the SRM which was proposed and analyzed in previous chapters for ordinary DAEs, especially for the index-2 DAEs (3.4), (3.5). Here we can take E = I even for 1 = 0 because (4.1) corresponds to (3.1) with B = GT . It is indicated in x3.1 that if we take 1 6= 0 then certain restrictions on choosing the initial iterate (cf. ( 3.14)) do not apply and, more importantly, the equation for xs is essentially not sti if the original problem is not sti . Hence, a non-sti time integrator can be used for any regularization parameter . For the Navier-Stokes application (4.2) we therefore choose 1 > 0 so that we can still take to be very small even when we use an explicit time discretization. So one SRM iteration is often good enough. However, we should not ignore the choice 1 = 0 because x3.1 also indicates that with this choice the computation can be particularly simple. For (4.2), when 1 > 0, although we use explicit time discretization, a symmetric positive de nite system relevant to the discretization of the operator I+ 1 graddiv still needs to be inverted. If we take 1 = 0, then we do not need to solve any system to obtain the discrete solution. In this case, (4.2) is not sti only for relatively large . So more than one SRM iterations are required generally. In the sequel, the convergence proof Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 85 and discretization stability analysis in x4.3 and x4.4 are mainly for the case of 1 > 0. The discussion for the case of 1 = 0 can essentially be carried out in a similar way. We will remark on this case in x4.3 and x4.4 and provide a numerical veri cation in x4.4. The importance of the treatment of the incompressibility constraint has long been recognized in the Navier-Stokes context. A classical approach is the projection method of [36], where one has to solve a Poisson equation for the pressure p with the zero Neumann boundary conditions which is, however, non-physical. Recently, a re-interpretation of the projection method in the context of the so-called pressure stabilization methods, or more generally, \pseudo compressibility methods" has been given in [97]. Some convergence estimates for the pressure can be obtained (cf. [101, 96]). In his review paper [97], Rannacher lists some well known examples of \pseudo-compressibility methods" (which are actually regularization methods): divu+ pt = 0; in [0; T ); pjt=0 = p0; (arti cial compressibility) divu+ p = 0; in [0; T ); (penalty method) divu p = 0; in [0; T ); @p @n j@ = p0 (pressure stabilization): If we generalize Baumgarte's stabilization to this PDAE example (4.1), we get ut + (u grad)u = u gradp+ f ; (4.3a) (divu)t + divu = 0: (4.3b) Eliminating ut from (4.3), we obtain an equation for p: p+ divu divf(u grad)u u fg = 0: We then nd that this stabilization can be seen as a kind of pressure stabilization with = 1. Although it works, since we do not have a singularity here, it still sets up a non-physical boundary condition for the Poisson equation for p. Also, in this formulation, equations for u and p are not uncoupled. Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 86 In the SRM formulation (4.2) we do not need to set up boundary conditions for p. So it should be more natural than various pressure-Poisson formulations. This method relates to the idea of penalty methods but, unlike them, the regularized problems are not sti for 1 > 0 or less sti for 1 = 0 since we can choose to be relatively large. Hence, more convenient (nonsti ) methods can be used for time integration, and nonlinear terms can be treated easily. We will indicate in x4.4 that has little to do with the stability of the discretization there, i.e. the stability restriction is satis ed for a wide range of for 1 > 0. We also indicate there that, in the case of small viscosity, the usual time step restrictions for the explicit schemes can be loosened. A similar procedure following [5] (Uzawa's iterative algorithm ) in the framework of optimization theory and economics has actually appeared in the Navier-Stokes context for the stationary Stokes equations (i.e. without the nonlinear term and the time-dependent term in (4.1)) with 1 = 0 using the augmented Lagrangian idea, see Fortin and Glowinski [50]. (Also see [54] for a related discussion.) Note that, in their procedure, 1 in (4.2c) is replaced by a parameter . They prove that = 1 is approximately optimal. For the nonlinear case, they combine Uzawa's algorithm with a linearization iteration. They claim convergence but nd it hard to analyze the convergence rate because their analysis depends on the spectrum of an operator which is non-symmetric in the nonlinear case. For the nonstationary case (4.1), the augmented Lagrangian method cannot be applied directly. Therefore, [50] rst discretizes (4.1a) with respect to the time t (an implicit scheme is used). Then the problem becomes a stationary one in each time step. Hence, Uzawa's algorithm can be applied and converges in each time step. So, for the nonstationary case, their iterative procedure is, in essence, to provide a method to solve the time-discretized problem. Thus, their iterative procedure has little to do with time-discretization, or in other words, they still do time-discretization directly for the problem (4.1). Consequently an implicit scheme is always appropriate because of the constraints Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 87 (4.1b), and a linearization is always needed to treat the nonlinear case. These properties are not shared by our method. We will prove that the convergence results of the SRM in the previous chapter still hold for the PDAE case (4.2). Hence, the solution sequence of (4.2) converges to the solution of (4.1) with the error estimate of O( s) after the sth iteration. Therefore, roughly speaking, the rate is about O( ). We prove the convergence results using the method of asymptotic expansions which is independent of the optimization theory and is also applicable to the steady-state case. In addition, when the nite element method is used, the di culty of constructing test functions in a divergence-free space to decouple the u; p system can be avoided by using the formulation of the SRM. We indicate here that, as many others do, we include the viscosity parameter in the error estimates. So the estimates could deteriorate when is very small. This is because we have an unresolved technical di culty, associated with our inability to obtain an appropriate upper bound for the nonlinear term and with the weaker elliptic operator u (which is dissipative) as ! 0. In the SRM formulation, a supplementary dissipative term 2graddivus is introduced without perturbing the solution. As indicated in [50] for the stationary case, the relative advantage of such methods may therefore become more apparent for small values of the viscosity. The chapter is organized as follows: In x4.2 we de ne some preliminaries and discuss regularity properties of the solution of (4.2). The convergence of the SRM for Navier-Stokes equations is proved in x4.3. Finally, in x4.4 a simple di erence scheme is discussed and some numerical experiments are presented. These numerical experiments are only exploratory in nature. To summarize, our objective in this chapter is to present a method for the nonstationary Navier-Stokes equations from the viewpoint of DAE regularization, and to provide a way to apply a DAE method to PDAEs. It appears that such a formulation is new in the Navier-Stokes context and it is worthwhile because: Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 88 Since need not be very small, the regularized problems in the sequence (4.2) are more stable/less sti . So more convenient di erence schemes, e.g. explicit schemes in time, can be used under theoretical assurance. If we take 1 > 0, this is also true for small . The problem of additional boundary conditions, which arises in the pressurePoisson formulation and projection methods, does not arise here. Finite element methods can be used easily and the elements do not have to conform to the incompressibility condition to separate the variables u and p. 4.2 Preliminaries and the Properties of the Regularized Problems Before we begin our analysis, we rst describe some notation and assumptions. As usual, we use Lp( ), or more simply Lp, to denote the space of functions which are pth-power integrable in , and kukp = (Z n Xi=1 upi dx) 1 p as its norm, where u = (u1; ; un). We denote the inner product in L2 by ( ; ) and let k k k k2. C1 is the space of functions continuously di erentiable any number of times in , and C10 consists of those members of C1 with compact support in . Hm is the completion in the norm kukHm = ( X 0 j j m kD uk2) 1 2 : We will consider the boundary conditions to be homogeneous, i.e. b 0 in (4.1c), to simplify the analysis. Nevertheless, through the inclusion of a general forcing term, the results may be generalized to the case of nonhomogeneous boundary values. We are interested in the case that (4.1) has a unique solution and the solution belongs to H2, where the arbitrary constant which the pressure p is up to is determined by Z p(x; ) dx = 0: (4.4) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 89 Hence, some basic compatibility condition is assumed (cf. [63]): aj@ = 0; diva 0: (4.5) Furthermore, we assume sup t2[0;T ]kfk M1; kakH2 M1; (4.6) where M1 is a positive constant. We take p0 in (4.2) satisfying (4.4). Hence, it is easy to see that ps satis es (4.4) for all s. For simplicity, we only consider the two-dimensional case. We can treat the threedimensional case in the same way, possibly with some more assumptions. Throughout the chapter M represents a generic constant which may depend on as we have explained in the introduction. We will also allow M to depend on the nite timeinterval T since we are not going to discuss very long time behavior in this chapter. At rst, we write down some inequalities: Poincar e inequality: kuk kgrad uk; if uj = 0: (4.7) More generally (see [87]), for u 2 H1( ) kuk C (kgrad uk+ j Z u dxj): (4.8) Young's inequality: abc 1 pap + 1 q bq + 1rcr (4.9) if a; b; c > 0; p; q; r > 1 and 1 p + 1q + 1r = 1. Holder's inequality: Z jf jjgjjhjdx kfkpkgkqkhkr (4.10) if p; q; r > 1 and 1 p + 1q + 1 r = 1. Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 90 Sobolev's inequality in two-dimensional space: kuk4 1 4 1 kuk 1 2kgrad uk 1 2 ; (4.11) where 1 = 2 if = R2. Suppose thatw stands for the di erence of two solutions of the SRM (4.2). Thenw satis es the homogeneous problem (4.27) (see next section). Hence, using the estimate in Lemma 4.2, uniqueness of the solution of the SRM (4.2) is easy to consider. The existence can be analyzed by following the standard existence argument for NavierStokes equations (e.g. [108, 62]) and that of penalized Navier-Stokes equations (e.g. [28]). Hence we assume the existence of the solution of the SRM and concentrate on the proof of the convergence of the method. Before we do so, we derive the following regularity results of the solution. Lemma 4.1 For the solution fus; psg of (4.2) there exists a constant 0 such that when 0 we have the following estimates kusk2H1 + Z T 0 ( 1 2 k(divus)tk2H1 + 2 2 kdivusk2H1 + k(us)tk2 + k usk2 + kpsk2H1) dt M [kak2H1 + Z T 0 (kfk2 + kgradps 1k2) dt]: (4.12) Proof: For simplicity of notation we denote us as v here. The proof for the case 1 = 0 is just the same as that in [28]. So we only consider the case 1 > 0. Hence, without loss of generality, we take 1 = 1 and 2 = . We then write (4.2) as vt 1 grad((divv)t + divv) + (v grad)v = v gradps 1 + f ; (4.13a) vj@ = 0;vjt=0 = a; (4.13b) ps = ps 1 1 ((divv)t + divv): (4.13c) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 91 The proof follows the ideas in [28]. Multiplying (4.13a) by v and integrating with respect to the space variables on the domain , we get 12 d dtkvk2 + 1 2 d dtkdivvk2 + kdivvk2 + kgradvk2 = ((v grad)v;v) (gradps 1;v) + (f ;v) 2 kdivvk2 + 1 2 kvk2kgradvk2 + 2 kvk2 + kgradps 1k2 + kfk2; where we use ((v grad)v;v) = 1 2((divv)v;v). Then let c = min( ; 2 ) and Y = kvk2 + 1 kdivvk2. Using Poincar e's inequality (4.7), we obtain d dtY + cY + 1 2( 1 Y )kgradvk2  (kfk2 + kgradps 1k2): (4.14) Note that Y (0) = kak2. Write (4.14) as d dt( 1 Y ) 1 2 kgradvk2( 1 Y ) 1 (kgradps 1k2 + kfk2): Applying a standard technique for solving linear di erential equations and taking appropriately small ( 0) so that 1 Y (0) 2 and 1 Z T 0 (kfk2 + kgradps 1k2) dt 4 ; we get 1 Y (t) 4 8t 2 [0; T ]: (4.15) Then, using the same technique and (4.14), we have Y kak2 exp( ct) +M exp( ct) Z t 0 (kfk2 + kgradps 1k2) exp(cz) dz M [kak2 + Z t 0 (kfk2 + kgradps 1k2) dz]: (4.16) Thus, (4.12) holds for kuk2. Integrating (4.14) directly and using (4.15) yields Z t 0 kgraduk2 dz M [kak2 + Z t 0 (kfk2 + kgradps 1k2) dz]: (4.17) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 92 To prove other estimates for (4.12) we de ne the operator Aw = 1 grad((divw)t + divw) w = g; (4.18) where w satis es wj@ = 0 and wjt=0 = a. Let q = 1 ((divw)t + divw): Then we have (noting divwjt=0 = 0) w + gradq = g; (4.19a) divw = Z t 0 q exp( (t z)) dz (4.19b) This is a general nonhomogeneous Stokes problem. Using the results described in [28] (or cf. [108]), we get k wk+ kgradqk M [kgk+ Z t 0 kgradqk exp( (t z)) dz]: (4.20) Applying Gronwall inequality, it is easy to obtain kgradqk M(kgk+ Z t 0 kgradqk dz) (4.21) and Z t 0 kgradqk2 dz M Z t 0 kgk2 dz =M Z t 0 kAwk2 dz: (4.22) It thus follows that k wk M(kgk+ Z t 0 kgk dz) =M(kAwk+ Z t 0 kAwk dz); (4.23) and then Z t 0 k wk2 dz M Z t 0 kgk2 dz =M Z t 0 kAwk2 dz: (4.24) From (4.19b) and (4.22), we thus have 1 2 Z t 0 kgraddivwk2 dz = Z t 0 kgradqk2 dz M Z t 0 kAwk2 dz: (4.25) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 93 Then 1 2 Z t 0 kgrad(divw)tk2 dz M Z t 0 kAwk2 dz (4.26) follows from (4.18). Now taking the scalar product of (4.13a) with Av, we have 1 2 k(divv)tk2 + 2 d dtkdivvk2 + 2 d dtkgradvk2 + kAvk2 = ((v grad)v; Av) (gradps 1; Av) + (f ; Av): Note that ((v grad)v; Av) kvk4kgradvk4kAvk 1 2 1 kvk 1 2kgradvkk vk 1 2kAvk (kAvk2+ k vk2) + 2 1 16 3 (kvk2kgradvk2)kgradvk2 [kAvk2+M2kAvk2 +M2 2(Z t 0 kAvk dz)2] + 2 1 16 3 (kvk2kgradvk2)kgradvk2; where we use (4.23) for the last inequality. Recall that we have estimates for kvk2 and R t 0 kgradvk2 dz. Therefore, taking (1 +M2) < 1 4, it is not di cult to obtain kgradvk+ 1 Z t 0 k(divv)tk2 dz + Z t 0 kAvk2 dz M [kgradak2 + Z t 0 (kfk2 + kgradps 1k2) dz + 2 Z t 0 kAvk2 dz] Taking such that M 2 < 1, we then get (4.12) for R t 0 kAvk2 dz and kgradvk. Noting (4.25), (4.26), from (4.2c), (4.12) also holds for R t 0 kgradpsk2 dz. Applying the inequality (4.8) and noting that ps satis es (4.4), yields the bound for R t 0 kpsk2 dz. Hence, from (4.2c) we can obtain (4.12) for R t 0 kdivvk2 dz and then R t 0 k(divv)tk2 dz. We thus complete the proof. 2 From this lemma, we see that if we choose p0 such that R t 0 kgradp0k2 dz is bounded then, by induction, all terms in the left of (4.12) are bounded for any given s. Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 94 4.3 Convergence of the SRM In this section, we estimate the error of the SRM (4.2) in the solution of (4.1) by using the asymptotic expansion technique as in the proof of Theorem 2.1 of Chapter 2 (see x2.8). Note that, in the Navier Stokes context, the asymptotic expansion method was used in [54] to obtain a more precise estimate for a penalty method for the stationary Stokes equations and in [108] to calculate a slightly compressible steady-state ow. We will mainly consider the case 1 > 0. Hence, we take 1 = 1 and 2 = for convenience. The result for 1 = 0 will be described in Remark 4.3. At rst we discuss a couple of linear auxiliary problems. Then we go to the proof. 4.3.1 Two linear auxiliary problems We discuss two linear problems in this section. One is wt grad(divw)t graddivw + (w grad)U + (V grad)w = w gradq + f ; (4.27a) wj@ = 0; wjt=0 = 0; (4.27b) where U, V and q are given functions. The other is wt + (V grad)w + (w grad)V = w gradp+ f; (4.28a) (divw)t + divw = g; (4.28b) wj@ = 0;wjt=0 = a; (4.28c) where V, g and a are given functions, a satis es the compatibility conditions (4.5) and g satis es (4.4). Now we show some properties of these two problems which will be used later in the proof of the convergence of SRM. Lemma 4.2 For the solution of problem (4.27), if U and V satisfy k k2H1 + Z T 0 k k2H2 dt M; (4.29) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 95 then we have the following estimate kwk2 + kdivwk2 M Z t 0 (kfk2 + kqk2) ds; (4.30a) kgradwk2 + Z t 0 ( kwtk2 + kdivwtk2) ds M Z t 0 (kfk2 + kqk2) ds: (4.30b) Proof: Multiplying (4.27a) by w and then integrating on the domain yields 1 2 d dtkwk2 + 1 2 d dtkdivwk2 + kdivwk2 + kgradwk2 = ((w grad)U;w) ((V grad)w;w) + (q; divw) + (f ;w) kgradUkkwk24 + 2kdivVkkwk24 + (q; divw) + (f ;w) (using (4.10)) 1 2 1 (kgradUk+ 1 2kdivVk)kwkkgradwk+ (q; divw) + (f ;w) (using (4.11)) 1 2 kgradwk2 + 1 2 (kgradUk+ 1 2kdivwk)2kwk2 + (q; divw) + (f ;w); where we have used ((V grad)w;w) = 2((divV)w;w). Therefore, we have d dt( kwk2 + kdivwk2) C(t)( kwk2 + kdivwk2) kgradwk2 ( + C(t))kdivwk2 + 2 (q; divw) + 2 (f ;w) kgradwk2 + 2kwk2 + 2 kfk2 + 1 kqk2 2 kfk2 + 1 kqk2; (4.31) where C(t) = 1 (kgradUk+ 1 2kdivVk)2: Noting that wjt=0 = 0 and divwjt=0 = 0, we thus get (4.30a). Now, multiplying (4.27) by wt, then integrating with respect to x over , we get kwtk2 + kdivwtk2 + 2 d dtkdivwk2 + d dtkgradwk2 = ((w grad)U;wt) + ((V grad)w;wt) + (q; divwt) + (f ;wt): (4.32) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 96 We use the inequalities listed in the previous section to estimate the right-hand side of (4.32) and have the following: ((w grad)U;wt) 4kwtk2 +M (kgradwk2 + kwk2); ((V grad)w;wt) 4kwtk2 +M sup jVj2kgradwk2; (f ;wt) 4kwtk2 +M kfk2; (q; divwt) 2kdivwtk2 +M kqk2; where the bounds of sup jUj and sup jVj can be obtained by using the inequality sup j j Mk k (see, e.g. [119]). Then, similarly to the procedure for obtaining (4.30a), we obtain (4.30b). 2 Next we consider problem (4.28). Lemma 4.3 There exists a solution for problem (4.28). Moreover, for the solution of (4.28), we have the following estimate: kwkH1 + Z T 0 (kwk2H2 + kwtk2 + kpk2H1) dt M (4.33) if R T 0 kfk2 dt and R T 0 kgk2H1 dt are bounded. Proof: First, we can solve divw from (4.28b) (noting that divwjt=0 = 0): divw = g1; (4.34) where g1 = exp( t) Z t 0 g exp( s)ds (4.35) satis es (4.4) since g does. By applying Corollary 2.4 in [54, p.23], the problem divw = g1; (4.36a) wj@ = 0 (4.36b) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 97 has many solutions. We pick one and denote it as wp. Then w := w wp satis es the linearized Navier-Stokes equations in the form of (4.28a) with a proper force term (denoted by f) and div w = 0; wj@ = 0 and wjt=0 = a wpjt=0: By noting that divwpjt=0 = g1jt=0 = 0, the basic compatibility conditions like (4.5) for w are satis ed. From (4.35) and the assumption for g, R T 0 (kg1k2H1+k(g1)tk2) dt is bounded. Hence, based on the estimates for the solution of (4.36) (see [1] and [54]), it is not di cult to get kwpkH1 + Z T 0 (k(wp)tk2 + kwpk2H2) dt M: (4.37) Thus R T 0 k fk2 dt is bounded. Simulating the regularity argument of [62] or [63] (multiplying the linearized Navier-Stokes equations by w, wt and P w where P is a projection operator (cf. [62]), respectively), we can obtain k wkH1 + Z T 0 (k wk2H2 + k wtk2 + kpkH1) dt M: (4.38) Therefore, (4.33) follows from (4.37) and (4.38). Using the estimate (4.38) and following a global existence argument (e.g. [62] or [108]), the existence of the solution for w can be obtained. We thus have the results of the lemma . 2 Remark 4.1 The uniqueness of the solution of (4.28) follows from the standard argument for Navier-Stokes equations (cf. [108]). 2 4.3.2 The error estimate of SRM In this section we prove the convergence of iteration (4.2) based on the same procedure described in the proof of Theorem 2.1 of Chapter 2 (see x2.8) . We describe our results in the following theorem. Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 98 Theorem 4.1 Let u and p be the solution of problem (4.1), and us and ps be the solution of problem (4.2) at the sth iteration. Then there exists a constant 0 such that when 0 we have the following error estimates for all t 2 [0; T ]: ku uskH1 M s; (4.39a) (Z T 0 kp psk2 dt) 1 2 M s; (4.39b) where T is any given nite number and s = 1; 2; . Proof: At rst, consider the case s = 1 of (4.2). Let u1 = u10 + u11 + + mu1m + Comparing the coe cients of like powers of , we thus have grad((divu10)t + divu10) = 0; (4.40a) grad((divu11)t + divu11) = (u10)t + (u10 grad)u10 u10 + gradp0 f ; (4.40b) grad((divu1i)t + divu1i) = (u1i 1)t + i 1 X j=1(u1j grad)u1i 1 j u1i 1 ; 2 i m+ 1; (4.40c) where (4.40a) satis es (4.2b) and (4.40b) and (4.40c) satisfy the homogeneous initial and boundary conditions corresponding to (4.2b). Now (4.40a) has in nitely many solutions in general. We should choose u10 not only to satisfy (4.40a) but also to ensure that the solution of (4.40b) exists. A choice of u10 is the exact solution u of (4.1), i.e. (u10)t + (u10 grad)u10 = u10 gradp + f ; (4.41a) (divu10)t + divu10 = 0; (4.41b) u10j@ = 0; u10jt=0 = a: (4.41c) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 99 Note that divu10jt=0 = diva = 0 and p is taken to satisfy (4.4). So u10 u and (4.40b) has the formgrad((divu11)t + divu11) = grad(p0 p): (4.42) Now we choose u11 and a corresponding p11 to satisfy (u11)t + (u10 grad)u11 + (u11 grad)u10 = u11 gradp11; (4.43a) (divu11)t + divu11 = p0 p; (4.43b) u11j@ = 0; u11jt=0 = 0: (4.43c) Again we have divu11jt=0 = 0 and let p11 satisfy (4.4). According to Lemma 4.3, u11 and p11 exist . Generally, supposing we have u1i 1, p1i 1 for i 2, choose u1i, p1i satisfying (u1i)t + (u10 grad)u1i + (u1i grad)u10 = u1i gradp1i i 1 X j=1(u1j grad)u1i 1 j ; (4.44a) (divu1i)t + divu1i = p1i 1; (4.44b) u1ij@ = 0; u1ijt=0 = 0; (4.44c) where we note that divu1ijt=0 = 0 and p1i satis es (4.4). Applying Lemma 4.3, all u1i and p1i; i = 0; 1; , exist and satisfy (4.33). Next we estimate the remainder of the asymptotic expansion after the (m+ 1)th power of . Denote u1m = u10 + u11 + + m+1u1m+1 (4.45) ( u1m also satis es (4.33)) and w1m = u1 u1m: (4.46) Then w1m satis es (w1m)t grad(divw1m)t graddivw1m Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 100 + (w1m grad)u1 + ( u1m grad)w1m = w1m m+2f(u1m+1)t + m+1 Xi=0 [(u1i grad)u1m+1 i] u1m+1g; (4.47a) w1mj@ = 0; w1mjt=0 = 0: (4.47b) Using regularity we have for u1i, u1m and u1 (see (4.12)) and Lemma 4.2 , we obtain kw1mk = O( m+1) and kgradw1mk = O( m+1). Therefore u1 = u10 + u11 + + mu1m +O( m+1): (4.48) in the H1-norm for the spatial variables. Noting u10 u, we thus obtain u1 u = O( ): (4.49) Furthermore, according to Lemma 4.2, we have kdivw1mk = O( m+ 32 ); (Z T 0 k(divw1m)tk2 dt) 1 2 = O( m+ 32 ): Then, by using (4.2c),(4.48),(4.41b), (4.43b), (4.44b) and the estimates for divw1m and (divw1m)t, it follows that p1 = p + p11 + + mp1m +O( m+ 12 ) (4.50) or p1 p = O( ): (4.51) in the sense of the L2-norm for both spatial and time variables, i.e. (R T 0 k k2 dt) 1 2 . Now we look at the second iteration s = 2 of (4.2). Let u2 = u20 + u21 + + mu2m + : Noting that (4.50) gives us a series expansion for p1 we obtain grad((divu20)t + divu20) = 0; (4.52a) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 101 grad((divu21)t + divu21) = (u20)t + (u20 grad)u20 u20 + gradp f ; (4.52b) grad((divu2i)t + divu2i) = (u2i 1)t + i 1 X j=1(u2j grad)u2i 1 u2i 1 gradp1i 1; 2 i m+ 1: (4.52c) Again, (4.52a) is combined with the initial and boundary conditions (4.2b), and (4.52b) and (4.52c) are combined with the corresponding homogeneous ones. As in the case of s = 1, we again choose u20 = u. We thus have grad((divu21)t + divu21) = 0: (4.53) Then u21 is constructed to satisfy (u21)t + (u20 grad)u21 + (u21 grad)u20 = u21 gradp21; (4.54a) (divu21)t + divu21 = 0; (4.54b) u21j@ = 0; u21jt=0 = 0: (4.54c) Obviously u21 = 0 and p21 = 0 is the solution of (4.54) and (4.4). In general, similarly to the case of s = 1, we choose u2i; p2i to satisfy (u2i)t + (u20 grad)u2i + (u2i grad)u20 = (4.55a) u2i gradp2i i 1 X j=1(u2j grad)u2i 1 j; (4.55b) (divu2i)t + divu2i = p1i 1 p2i 1; (4.55c) u2ij@ = 0; u2ijt=0 = 0 (4.55d) for 2 i m + 1, where p2i satis es (4.4). By the same procedure as for s = 1 we obtain error equations similar to (4.47) with the addition of a remainder term grad(p1 p1m) in the right-hand side, where p1m stands for the asymptotic expansion (4.50) of p1. Applying Lemma 4.2 again, we get u2 = u20 + u21 + + mu2m +O( m+1): (4.56) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 102 Noting u21 0, u2 u = O( 2): (4.57) Then, using (4.2b),(4.56),(4.54b) and (4.55c), we conclude p2 = p + p21 + + mp2m +O( m+ 12 ) (4.58) or p2 p = O( 2): (4.59) since p21 0. We can repeat this procedure, and, by induction for s (choosing m larger than s), conclude the results of the theorem. 2 Remark 4.2 Corresponding to Theorem 2.1, we expect that the error estimates (4.39) also hold for the SRM (4.2) with 1 = 0, at least, away from t = 0. 2 Remark 4.3 In Theorem 4.1, we nd that the result for p is in a weaker norm R T 0 k k2 dt. This is because we have di culty in estimating the rst order timederivative of the right-hand side of (4.47), or concretely, the term R T 0 k(u1m+1)ttk2 ds. In [63] (Corollary 2.1) it is shown that R T 0 k(u1m+1)ttk2 ds may be unbounded as t! 0 if we only assume the local compatibility conditions (4.5). In the case that this integral is bounded for 0 < t < T , we can get kp psk+ (Z t 0 k(p ps)tk2 ds) 12 M s: (4.60) Otherwise, we only can expect that (4.60) holds away from t = 0 by following the argument in [63]. 2 Remark 4.4 Multiplying (4.27) by Aw, where A is the operator de ned by (4.18), and following the later steps of the proof of Lemma 4.1, we can get Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 103 1 2 Z T 0 (kgrad(divw)tk2 + kgraddivwk2) dt M Z T 0 (kfk2 + kgrad qk2) dt: Using this result to estimate the remainders of the asymptotic solutions in the proof of Theorem 4.1, we obtain (Z T 0 kp psk2H1 dt) 1 2 M s: (4.61) 24.4 Discretization Issues and Numerical Experiments In previous sections, we have proposed the SRM and performed some basic analysis on it. The SRM yields a sequence of PDEs which are to be solved numerically. The problem at the sth iteration can be written as: (us)t grad( 1(divus)t + 2divus) + (us grad)us = us + rs; (4.62a) usj@ = 0;usjt=0 = a; (4.62b) where rs(t) is the known inhomogeneity rs = gradps 1 + f : (4.63) A variational formulation of (4.62) gives: Find us 2 H10 such that d dt(us; ) + 1 d dt(divus; div ) + 2(divus; div ) + (gradus;grad ) + b(us;us; ) = (rs; ); 8 2 H10; (4.64a) usjt=0 = a ; divusjt=0 = 0; (4.64b) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 104 where the trilinear form b(u;v;w) = ((u grad)v;w): From (4.64) we see that the nite element method in the spatial variables, combined with time discretizations, can be easily adopted. Note that we do not need to construct divergence-free test functions to separate the variables u and p. Nevertheless, in this section, we are not going to discuss nite element methods further. Some discussions on using the SRM with the nite element method will be given in the next chapter for a problem in reservoir simulations. Some numerical experiments are also given there. Here we only consider a very simple rst-order di erence scheme (forward Euler scheme in the time direction) in two dimensional space, as an initial attempt towards the discretization of the sequential regularization method for the PDAE. Concretely, we consider a rectangular domain such that an equidistant mesh can be used. Let (u; v)T stand for the approximation of us, and let k; hx; hy denote step sizes in time and spatial direction, respectively. Without loss of generality, we assume that hx = hy = h and that the domain is a unit square. Thus, mesh points can be expressed as xi = ih; i = 0; 1; ; I; yj = jh; j = 0; 1; ; ; J ; tn = nk; n = 0; 1; ; N; N = [T=k]: The di erence scheme reads: u _ t 1(u x _ x + v y _ x) _ t = 2(u x _ x + u y _ x) (uu x + vu y) + (u x _ x + u y _ y) + ru; (4.65a) v _ t 1(u x _ y + u y _ y) _ t = 2(u x _ y + u y _ y) (uv x+ vv y) + (v x _ x + v y _ y) + rv; (4.65b) uj@ = 0; vj@ = 0; ujt=0 = au; vjt=0 = av; (4.65c) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 105 where u = uni;j; u _ t = un+1 i;j uni;j k ; u _ x = uni+1;j uni;j h ; u x = uni;j uni 1;j h : u _ y and u y can be de ned accordingly and the de nitions for v are similar. Obviously, this is a rst-order scheme explicit in time, where the nonlinear term is discretized somewhat arbitrarily. The scheme is easy to implement. Next we discuss its stability. For simplicity, we analyze the linear case (corresponding to the Stokes equations) rst, and consider the full nonlinear equations (4.65) in Remark 4.7 below. We write the linear case of (4.65) as follows : u _ t 1(u x _ x + v y _ x) _ t = 2(u x _ x + u y _ x) + (u x _ x + u y _ y) + ru; (4.66a) v _ t 1(u x _ y + u y _ y) _ t = 2(u x _ y + u y _ y) + (v x _ x + v y _ y) + rv; (4.66b) uj@ = 0; vj@ = 0; ujt=0 = au; vjt=0 = av: (4.66c) Here we take 1 = 1 and 2 = . The result for the case of 1 = 0 will be given in Remark 4.6. The following theorem gives the stability estimate for (4.66) in the sense of the discrete L2-norm: kwhk2h = h2 I 1 Xi=0 J 1 X j=0(wi;j)2; (4.67) where wh = (wi;j); i = 0; 1; ; I 1; j = 0; 1; ; J 1. Theorem 4.2 Let u and v be the solution of (4.66) and A = (kuk2h + kvk2h) + ku x + v yk2h + (ku xk2h + ku yk2h + kv xk2h + kv yk2h): (4.68) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 106 If k h2 1 c, where c is any constant in (0; 1), then A+ k N 1 X n=0 k(u x + v y) _ tk2h M max 0 tn T(kruk2h + krvk2h) (4.69) where M is a generic constant dependent on and c. Proof: We rst write down the following identities and inequalities: some di erence identities [75]: ( ) x = x + xE 1 x ; (4.70a) ( ) _ x = _ x + _ xE1 x ; (4.70b) 2 _ t = ( 2) _ t k( _ t)2; (4.70c) _ x x = ( _ x) x ( x)2; (4.70d) where the translation operator Ei x (x; y; t) = (x+ ih; y; t). an di erence inequality [75]: hk xkh 2k kh (4.71) a discrete version of the Poincar e inequality (cf. [66]): k k2h k xk2h + k yk2h (4.72) if satis es homogeneous boundary conditions. Multiplying (4.66a) by au+bu _ t and (4.66b) by av+bv _ t and adding, then summing for all (i; j); i = 1; ; I 1; j = 1; ; J 1 where we use the di erence identities (4.70) (omitting lots of tedious algebraic manipulations), we obtain Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 107 a(kuk2h + kvk2h) _ t + 1 2( b+ a)(ku x + v yk2h) _ t + 1 2 b(ku xk2h + ku yk2h + kv xk2h + kv yk2h) _ t + (b ak)(ku _ tk2h + kv _ tk2h) + a ku x + v yk2h + (b 1 2(b + a)k)k(u x + v y) _ tk2h + a(ku xk2h + ku yk2h + kv xk2h + kv yk2h) 1 2 bk(ku x _ tk2h + ku y _ tk2h + kv x _ tk2h + kv y _ tk2h) M (kruk2h + krvk2h) + 1 2 a(kuk2h + kvk2h) + 1 2 b (ku _ tk2h + kv _ tk2h); where > 0 can be chosen less than c=b. Applying (4.71) and (4.72), we get h2((ku x _ tk2h + ku y _ tk2h + kv x _ tk2h + kv y _ tk2h) 4(ku _ tk2h + kv _ tk2h) and kuk2h + kvk2h ku xk2h + ku yk2h + kv xk2h + kv yk2h; respectively. Then we can choose a and b such that b ak k h2 12b > 0; b 1 2(b + a)k > 0 and obtain (A) _ t + d A+ k(u x + v y) _ tk2h M (kruk2h + krvk2h); where d is a constant independent of k, h, and . From this inequality, it is not di cult to see that (4.69) holds. 2 Remark 4.5 From (4.69) of Theorem 4.2, we nd that the value of will not a ect the stability of the di erence scheme. This means that the forward Euler scheme in the time direction works for any value of . Also, the time step restriction k (1 c)h2= is actually loosened in the case of small viscosity (or large Reynolds number) which people are often interested in. This implies that the explicit scheme (4.66) to which an appropriate discretization of the nonlinear term (see the next remark) is added works very well. It enables us not only to avoid the complicated iteration procedure for nonlinear equations but also to choose the time step fairly widely in the case of small viscosity. 2 Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 108 Remark 4.6 We have mentioned before that sometimes we may like to take 1 = 0 to avoid solving any algebraic system. Following the same procedure as the proof of Theorem 4.2, we can get the stability condition for the case of 1 = 0. That is, k m h2, where m is a positive constant independent of ,h and . We thus see that the stability of (4.66) with 1 = 0 depends on the parameter . This coincides with our experience with sti problems discretized by explicit schemes. Fortunately, using the SRM, we do not need to take very small. So the time step restriction is not much worse than the usual one corresponding to an explicit scheme applied to a non-sti problem. 2 Remark 4.7 For the nonlinear case (4.65), when the viscosity is not small, we expect similar results since the nonlinear term can be dominated by the viscous term. When the viscosity is small, however, the scheme (4.65) is unstable. Although numerical computations indicate that we do get better stability if we increase 2, i.e. some kind of dissipation e ect is obtained (we must note that such dissipation becomes small when the incompressibility condition is close to being satis ed), we suggest using spatial discretizations with better stability properties, e.g. upwinding schemes (cf. [102]), in the case of small viscosity. 2 Remark 4.8 Applying corresponding di erence identities for a nonuniform mesh (see e.g. [106]), the results of Theorem 4.2 may be generalized to di erence schemes (4.66) on a nonuniform mesh. Hence, the di erence scheme may be used for problems de ned on more general domains. 2 Next we explain our theoretical results by calculating the solution of an arti cial example. Example 4.1 Consider the Navier-Stokes equations (4.1) with the exact solution Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 109 u = (u; v)T u = 50x2(1 x)2y(1 y)(1 2y)[1 + exp( t)]; v = 50y2(1 y)2x(1 x)(1 2x)[1 + exp( t)]; p = [ x(x2 + 2) y(y2 2) + 1 3][1 + exp( t)]2 As indicated in x2.6 of Chapter 2, to carry out the SRM iterations, we do not need to store the entire approximation of ps 1 on [0; T ] for calculating us. Assuming that the number of the SRM iterations is chosen in advance, we can rearrange the computational order to make the storage requirements independent of N , where N represents the number of the mesh lines in the t direction. We rst use constant steps k = 0:01 and h = 0:1. At a given time t, we use `eu0 to denote the absolute discrete L2-error in us while `ep0 denotes the absolute discrete L2-error in ps. Table 4.1 summarizes the computational results of the di erence scheme (4.65) with 1 = 2 = 1 and viscosity = 0:1. iteration error at ! t = k t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 5e-1 1 eu 4.65e-3 2.69e-1 1.55e-1 1.31e-1 1.15e-1 1.08e-1 ep 2.49e-1 1.96e-1 1.57e-1 1.36e-1 1.25e-1 1.21e-1 2 eu 2.16e-3 2.53e-2 3.12e-2 3.18e-2 3.12e-2 3.06e-2 ep 1.80e-1 9.28e-2 7.37e-2 6.74e-2 6.48e-2 6.35e-2 3 eu 2.15e-3 1.77e-2 2.28e-2 2.48e-2 2.55e-2 2.57e-2 ep 1.80e-1 8.81e-2 6.69e-2 6.10e-2 5.91e-2 5.83e-2 1e-3 1 eu 2.14e-3 1.73e-2 2.21e-2 2.41e-2 2.48e-2 2.50e-2 ep 1.80e-1 8.78e-2 6.61e-2 6.01e-2 5.82e-2 5.75e-2 Table 4.1: SRM errors for = 0:1 without upwinding We notice that the errors improve as the iteration proceeds until s reaches the discretization accuracy O(h), where s is the number of iterations. For small viscosity, say = 0:001, the di erence scheme (4.65) does not work. The errors blow up around t = 1. When we increase 2, say to 50, we do get pretty good results around t = 1; however, the errors still blow up at a later time. This Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 110 suggests that the scheme is not stable for small viscosity . So we next discretize the nonlinear term using the upwinding scheme given in [102]. For the case of small viscosity, e.g. = 0:001, we get good results (see Table 4.2). iteration error at ! t = k t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 5e-1 1 eu 4.66e-3 2.26e-1 2.50e-1 2.29e-1 2.13e-1 2.03e-1 ep 2.58e-1 1.06e-1 6.67e-2 5.45e-2 5.12e-2 5.04e-2 2 eu 2.16e-3 7.74e-2 8.78e-2 9.13e-2 9.34e-2 9.53e-2 ep 1.84e-1 8.81e-2 6.22e-2 5.39e-2 5.11e-2 5.02e-2 3 eu 2.14e-3 7.69e-2 8.71e-2 9.06e-2 9.29e-2 9.48e-2 ep 1.83e-1 8.78e-2 6.21e-2 5.39e-2 5.11e-2 5.01e-2 1e-3 1 eu 2.14e-3 7.69e-2 8.72e-2 9.07e-2 9.29e-2 9.49e-2 ep 1.83e-1 8.78e-2 6.21e-2 5.39e-2 5.11e-2 5.01e-2 Table 4.2: SRM errors for = 0:001 with upwinding Recall that according to Remark 4.5, in the case of small viscosity, the time step size can be increased to some extent without adverse stability e ects. To demonstrate this, we take k = h = 0:1, and = 0:001. The numerical results in Table 4.3 support our claim. iteration error at ! t = k t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 1e-3 1 eu 2.18e-2 8.61e-2 9.43e-2 9.70e-2 9.86e-2 9.99e-2 ep 1.83e-1 8.83e-2 6.26e-2 5.42e-2 5.13e-2 5.03e-2 Table 4.3: SRM errors for = 0:001 with a pretty large time step k = h = 0:1 Although we use explicit schemes for SRM (4.2) with 1 > 0, we still have to solve a banded symmetric positive de nite system. An alternative is to take 1 = 0 to avoid solving any algebraic systems. Table 4.4 shows the computational results of the di erence scheme (4.65) with 1 = 0 and 2 = 1. We take viscosity = 0:1, h = 0:1 and k = 0:0005. Good results are obtained except for the pressure near t = 0 (cf. Remark 4.4) . Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 111 iteration error at ! t = k t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 5e-1 2 eu 5.64e-3 3.57e-2 2.94e-2 2.71e-2 2.62e-2 2.60e-2 ep 2.92e-0 9.70e-2 7.03e-2 6.17e-2 5.87e-2 5.77e-2 Table 4.4: SRM errors for = 0:1 with 1 = 0 Chapter 5 SRM for the Simulation of Miscible Displacement in Porous Media 5.